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Abstract. The conformal formulation of the Einstein constraint equations is 
first reviewed, and wc then consider the design, analysis, and implementation 
of adaptive multilevel finite elemcnt-tyijc immcrical methods for the resulting 
coupled nonlinear elliptic system. VVc derive weak formulations of the coupled 
» ' I constraints, and review some new developments in the solution theory for the 

(2J[^ constraints in the cases of constant mean extrinsic curvature (CMC) data, near- 

I i CMC data, and arbitrarily prescribed mean extrinsic curvature data. We then 

outline some recent results on a priori and a posteriori error estimates for a broad 
CO class of Galerkin-type approximation methods for this system which includes 

^ techniques such as finite element, wavelet, and spectral methods. We then 

use these estimates to construct an adaptive finite element method (AFEM) for 
solving this system numerically, and outline some new convergence and optimality 
y—{ results. We then describe in some detail an implementation of the methods 

using the FETK software package, which is an adaptive multilevel finite element 
• code designed to solve nonlinear elliptic and parabolic systems on Ricmannian 

'"^ manifolds. We finish by describing a simplex mesh generation algorithm for 

compact binary objects, and then look at a detailed example showing the use 
00 of FETK for numerical solution of the constraints. 
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1. Introduction 

One of the primary goals of numerical relativity is to compute solutions of the full 
Einstein equations and to compare the results of such calculations to the expected data 
from the current or next generation of interferometric gravitational wave observatories 
such as LIGO. Such observatories are predicted to make observations primarily of 
"burst" events, mainly supernovae and collisions between compact objects, e.g., black 
holes and neutron stars. Computer simulation of such events in three space dimensions 
is a challenging task; a sequence of reviews that give a fairly complete overview of 
the state of the field of numerical relativity at the time they appeared are [i , 3]. 
While there is currently tremendous activity in this area of computational science, 
even the mathematical properties of the Einstein system are still not completely 
understood; in fact, even the question of which is the most appropriate mathematical 
formulation of the constrained evolution system for purposes of accurate longtime 
numerical simulation is still being hotly debated (cf. [4, •)]). 

As is well known, solutions to the Einstein equations are constrained in a manner 
similar to Maxwell's equations, in that the initial data for a particular spacetime 
must satisfy a set of purely spatial equations which are then preserved throughout 
the evolution (cf. [ ]). As in electromagnetism, these constraint equations may be 
put in the form of a formally elliptic system; however for the Einstein equations, the 
resulting system is a set of four coupled nonlinear equations which are generally non- 
trivial to solve numerically; moreover, the solution theory is only partially understood 
at the moment (cf. [7, 8, 9]). Until recently, most results were developed only in 
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the case of constant mean extrinsic curvature (CMC) data, leading to a decoupling 
of the constraints that allows them to be analyzed independently. A complete 
characterization of CMC solutions on closed manifolds appears in [10]; see also the 
survey [1 1] and the recent work on weak solutions [12, l.'J]. Of the very few non-CMC 
results available for the coupled system are those for closed manifolds in [7] ; these and 
all other known non-CMC results [14, 15] hold only in the case that the mean extrinsic 
curvature is nearly constant (referred to as near-CMC). Some very recent results on 
weak and strong solutions to the constaints in the setting of both CMC and near-CMC 
solutions on compact manifolds with boundary appear in [16]. Related results on weak 
and strong solutions to the constaints in the setting of CMC, near-CMC, and truly 
non-CMC (far-from-CMC) solutions on closed manifolds appear in [17, 18]. 

Numerous efforts to develop effective numerical techniques for the constraint 
equations, and corresponding high-performance implementations, have been 
undertaken over the last twenty years; the previously mentioned numerical relativity 
reviews [1, 2, .3] give a combined overview of the different discretization and solver 
technology developed to date. A recent review that focuses entirely on the constraints 
is [19]; this work reviews the conformal decomposition technique and its more recent 
incarnations, and also presents an overview of the current state of numerical techniques 
for the constraints in one of the conformal forms. While most previous work has 
involved finite difference and spectral techniques, both non-adaptive and adaptive, 
there have been previous applications of finite element (including adaptive) techniques 
to the scalar Hamiltonian constraint; these include [20, 21, 22, 23, 24]. An initial 
approach to the coupled system using adaptive finite element methods appears in [.. ]. 
A complete theoretical analysis of adaptive finite element methods (AFEM) for a 
general class of geometric PDE appears in [2G, 27], and a complete analysis of AFEM 
for the Einstein constraints appears in [-' ], including proofs of convergence and 
optimality. 

While yielding many useful numerical results and new insights into the Einstein 
equations, most of the approaches previously used for the constraints suffer from one 
or more of the following limitations: 

• The 3-metric must be conformally flat; 

• The extrinsic curvature tensor must be traceless or vanish altogether; 

• The momentum constraint must have an exact solution; 

• The problem must be spherically symmetric or axisymmctric; 

• The domain must be covered by a single coordinate system which may contain 
singularities; 

• The conformal metric must have a scalar curvature which can be easily computed 
analytically; 

• The gauge conditions and/or source terms must be chosen so that the constraints 
decouple, become linear, or both; 

• There is very little control of the approximation error, or even a rigorously derived 
a priori error estimate; 

• One must have access to a large parallel computer in order to obtain reasonably 
accurate solutions on domains of physical interest. 

In this paper, we describe a general approach based on adaptive finite element methods 
that allows one to avoid most of the limitations above. 
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More precisely, we develop a class of adaptive finite element methods for producing 
high-quality numerical approximations to solutions of the constraint equations in 
the setting of general domain topologies and general non-constant mean curvature 
coupling of the two constraint equations. Our approach is based on the adaptive 
approximation framework developed in [26, 27, 28], which is a theoretical framework 
and corresponding implementation for adaptive multilevel finite element methods 
for producing high-quality reliable approximations to solutions of nonlinear elliptic 
systems on manifolds. As in earlier work, we employ a 3-t-l splitting of spacetime and 
use the York conformal decomposition formalism to produce a covariant nonlinear 
elliptic system on a 3-manifold. This elliptic system is then written in weak form, 
which offers various advantages for developing solution theory, approximation theory, 
and numerical methods. In particular, this allows us to establish a priori error 
estimates for a broad class of Galerkin-type methods which include not only finite 
element methods, but also wavelet-based methods as well as spectral methods. Such 
estimates provide a base approximation theory for establishing convergence of the 
underlying discretization approach. We also derive a posteriori error estimates which 
can be used to drive adaptive techniques such as local mesh refinement. We outUne 
a class of simplex-based adaptive algorithms for weakly-formulated nonlinear elliptic 
PDE systems, involving error estimate-driven local mesh refinement. 

We then describe in some detail a particular implementation of the adaptive 
techniques in the software package FETK [26], which is an adaptive multilevel finite 
element code based on simplex elements. This software is designed to produce provably 
accurate numerical solutions to a large class of nonlinear covariant elliptic systems of 
tensor equations on 2- and 3-manifolds in an optimal or nearly-optimal way. It employs 
a posteriori error estimation, adaptive simplex subdivision, unstructured algebraic 
multilevel methods, global inexact Newton methods, and numerical continuation 
methods for the highly accurate numerical solution of nonlinear covariant elliptic 
systems on 2- and 3-manifolds. The FETK implementation has several unusual 
features (described in Section 6) which make it ideally suited for solving problems 
such as the Einstein constraint equations in an adaptive way. Applications of FETK 
to problems in other areas such as biology and elasticity can be found in [29, 30, 31]. 

Outline of the paper 

We review the classical York conformal decomposition in Section 2.2. In Section 2.3 
we give a basic framework for deriving weak formulations. In Section 2.4 we briefly 
outline the notation used for the relevant function spaces. In Section 2.5 we go over a 
simple weak formulation example. In Section 2.6 we derive an appropriate symmetric 
weak formulation of the coupled constraint equations, and summarize a number of 
basic theoretical results. In Section 2.7 we also derive the linearized bilinear form 
of the nonlinear weak form for use with stability analysis or Newton-like numerical 
methods. A brief introduction to finite element methods for nonlinear elliptic systems 
is presented in Section 3.1. Adaptive methods are described in Section 3.3, and 
residual-type error indicators are derived in Section 3.4. A derivation of the a 
posteriori error indicator for the constraints is give in Section 3.5. We give two 
give an overview of a priori error estimates from [26, 27, 28] for general Galerkin 
approximations to solutions equations such as the momentum and Hamiltonian 
constraints in Section 3.2. The numerical methods employed by FETK are described 
in detail in Section 6, including the finite element discretization, the residual-based 
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a posteriori error estimator, the adaptive simplex bisection strategy, the algebraic 
multilevel solver, and the Newton-based continuation procedure for the solution of 
the nonlinear algebraic equations which arise. Section 6 describes a mesh generation 
algorithm for modeling compact binary objects, outlines an algorithm for computing 
conformal Killing vectors, describes the numerical approximation of the ADM mass, 
and gives an example showing the use of FETK for solution of the coupled constraints 
in the setting of a binary compact object collision. The results are summarized in 
Section 7. 



2. The Constraint Equations in General Relativity 

2.1. Notation 

Let (A1,7ab) be a connected compact Riemannian d-manifold with boundary 
(Tab), where the boundary metric dob is inherited from ^ab- In this paper we are 
interested only in the Riemannian case, where we assume ^ab is strictly positive a.e. in 
M.. Later we will require some additional assumptions on ^ab such as its smoothness 
class. To allow for general boundary conditions, we will view the boundary (c? — 1)- 
submanifold dJ^ (which we assume to be oriented) as being formed from two disjoint 
submanifolds OqM and diM, i.e., 

d^M^JdiM^dM, ^nMC^^lM = %. (1) 

When convenient in the discussions below, one of the two submanifolds d^M or diM 
may be allowed to shrink to zero measure, leaving the other to cover dAi. An 
additional technical assumption at times will be non-intersection of the closures of 
the boundary sets: 



doM n diM = 0, (2) 

This condition is trivially satisfied if either d^M or diM shrinks to zero measure. It is 
also satisfied in practical situations such as black hole models, where d^AA represents 
the outer-boundary of a truncated unbounded manifold, and where diM represents 
the inner-boundary at one or more black holes. In any event, in what follows it 
will usually be necessary to make some minimal smoothness assumptions about the 
entire boundary submanifold dM , such as Lipschitz continuity (for a precise definition 
see [;-!2]). 

We will employ the abstract index notation (cf. [ ]) and summation convention 
for tensor expressions below, with indices running from 1 to d unless otherwise noted. 
Covariant partial differentiation of a tensor t""^ "''bi - b using the connection provided 

by the metric 7ab will be denoted as t""^ "''bi - b c '^'^ Dct"^^ °''bi- -6 ■ Denoting the 
outward unit normal to dAi as rib, recall the IDivergence Theorem for a vector field 
on M (cf. [:-! ']): 

w^.f, dx — w^Ub ds, (3) 
/M ' JdM 

where dx denotes the measure on A4 generated by the volume element of 7ab: 



dx = \/det "fab dx^ ■ ■ ■ dx'^, (4) 

and where ds denotes the boundary measure on dJ^ generated by the boundary 
volume element of aab. Making the choice w'' = Uai...afcW°i - °'=^ and forming the 
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divergence 'u/'.i^ by applying the product rule leads to a useful integration-by-parts 
formula for certain contractions of tensors: 

JM JdM Jm 

When fc = this reduces to the familiar case where u and v are scalars. 
2.2. The York Decomposition 

As discussed in the introduction, the most common form of the initial data 
problem for numerical work is the coupled elliptic system obtained from the York 
decomposition [35, 36] . We employ the standard notation whereby the spatial 3-metric 
is denoted jab with indices running from 1 to 3. The Hamiltonian and momentum 
constraints are 

R + {trKf - KabK"'' = 167rp (6) 

and 

DbiK"'' ~ -f'^hrK) 8nf, (7) 

where R is the scalar curvature of jabi is the extrinsic curvature of the initial 
hypersurface, and p and j'" are the mass density and current. 

As is well known, (6) and (7) form an under-determined system for the 
components of the 3~metric and extrinsic curvature. The York conformal 
decomposition method identifies some of these components as "freely specifiable," i.e., 
source terms similar to the matter terms p and j", and the remaining components as 
constrained. The main strength of the formalism is that it does this in such a way that 
the constraint system becomes a coupled set of nonlinear elliptic equations in these 
components. In only the most general setting do the fully coupled equations have to 
be solved; in many applications they can be decoupled and further simplification of the 
problem may eliminate one or more of the non-linearities. The primary disadvantage 
of the formalism is the difficulty of controlling the physics of the initial data generated; 
the matter terms as well as the unconstrained parts of and Kab are not related 
directly to important physical and geometrical properties of the initial data. One 
therefore generally relies on an evolution, i.e., a computation of the entire spacetime, 
to determine "what was in" the initial data to begin with. Some of the difficulties 
with use of the formalism have been recently overcome; cf. [19] for a recent survey. 

The formalism is summarized in a number of places (see, e.g., [35, 30, Ifl]), in this 
section we give only a brief description of the most standard form. The manifold A4 
is endowed with a Riemannian metric 7ah and the solution to the initial data problem 
is, in part, a metric jab related to jab by 

lab = i'^lab- (8) 

In what follows all hatted quantities are formed out of jab in the usual way, e.g., the 
covariant derivative Da, the Riemann tensor Rabcd, etc., while unhatted quantities are 
formed out of jab- 

There are four constraint equations and 12 components of jab and Kab. The 
eight freely specifiable components consist of the conformal part of the 3-metric and 
the trace and the trans verse-traceless parts of the extrinsic curvature. The remaining 
constrained parts are the conformal factor, 0, and the "longitudinal potential," W^, 
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of the extrinsic curvature. Specifically, we decompose the extrinsic curvature tensor 
as 

j^ab ^ + {LWy^) + -(jy-^'^hrK, (9) 

3 

3 

(LVF)"'' = D^'W^ + D^W - '^T^DcW", (12) 

3 

where trif = Y'-'^Kab- Note that A"-'' and (LW)'^^ are traceless by construction and 

is a freely specifiable transverse-traceless tensor. The matter terms are decomposed 

as 

P^pr', r = j"r'°. (13) 

In these variables the Hamiltonian constraint (6) becomes 

= + ^(trif)205 - + {LW)abf(t>'' - 271/50-3, (14) 

where Acj) — DaD°'4' ^^^d where {TabY = T"'''Tab, following the notation introduced 
in (34). In these variables the momentum constraint (7) becomes 

DbiLW)'''' = ^/l)"trif + Sirf. (15) 

In this article, we are mainly interested in formulations of the constraints on 
manifolds with boundary. The primary motivation for this is the desire to develop 
techniques for producing high-quality numerical solutions to the constraints, which 
require the use of mathematical formulations of the constraints involving finite 
domains with boundary. In order to completely specify the strong (and later, the 
weak) forms of the constraints on manifolds with boundary, we need to specify the 
boundary conditions. This is very problem dependent; however, we would like to at 
least include the case of the vector Robin condition for asymptotically flat initial data 
given in [''"]. This is 

{LW)''^h, (^<5? - ^-h^h,^ + ^W" (^5? - ^-h^h,^ = 0{R-') (16) 

where R is the radius of a large, spherical domain. The right hand side could be taken 
to be zero. Noting that 5^ + n°'ni, is the inverse of — 1/2 h"-hb we may compute 

{LWr'h, + Aiy" (^5- + 3 ^ p ^^^^ 

Hence we will consider the following linear Robin-like condition, which is general 
enough to include (17) and more recently proposed boundary conditions [ ]: 

{LWy^fib + C^W^ = Z'' ondiM. (18) 

Similarly, we are interested in analyzing the case of a Robin-like boundary condition 
with the Hamiltonian constraint: 

hab''(t> + c(j)= z ondiM. (19) 

Equations (14)-(15) are known to be well-posed only for restricted problem data 
and manifold topologies [39, 10, 40, 41, 42, 37, 43, 7, 44, 8, 9]; most of the existing 
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results are for the case of constant mean extrinsic curvature (CMC) data on closed 
(compact without boundary) manifolds, with some results for near-CMC data [ ]. 
Some very recent results on weak and strong solutions to the constaints in the setting 
of CMC and near-CMC solutions on compact manifolds with boundary appear in 
Related results on weak and strong solutions to the constaints in the setting of CMC, 
near-CMC, and truly non-CMC (far from-CMC) solutions on closed manifolds appear 
in [17, 18]. 



2.3. General Weak Formulations of Nonlinear Elliptic Systems 

Consider now a general second-order elliptic system of tensor equations in strong 
divergence form over a Riemannian manifold M. with boundary: 



A' 



\x\u\u\^).,a + B\x\u^ 



u\x'') 



mM, (20) 
on diM, (21) 
E\x^) ondoM, (22) 



where 



1 < i, j,k < n, 

nd ^ ^nd^ B -.Mx 



r>nd 



(23) 
(24) 
(25) 



1 < a,b,c < d, 
A:MxWx 

C : diM X M" X E : doM x 

The divergence-form system (20)-(22), together with the boundary conditions, can be 
viewed as an operator equation 

F(u) = 0, E:Bi<^B*2, (26) 

for some Banach spaces Bi and B2, where B2 denotes the dual space of B2. 

Our interest here is primarily in coupled systems of one or more scalar field 
equations and one or more d- vector field equations. The unknown n- vector then in 
general consists of scalars and d-vectors, so that n = + ■ d. To allow the 
n-component system (20)-(22) to be treated notationally as if it were a single n-vector 
equation, it will be convenient to introduce the following notation for the unknown 
vector and for the metric of the product space of scalar and vector components of 



(27) 





" 7ll^ 
























("=) 

Tab 







If 



is a d- vector we take 7^^^ — jab] if is a scalar we take j'^l' — 1 



(fc) 



The weak form of (20)-(22) is obtained by taking the L^-inner-product between 
a vector (vanishing on doM) lying in a product space of scalars and tensors, and 
the residual of the tensor system (20), yielding: 



dx = 0. 



(28) 



M 



Due to the definition of Gij in (27), this is simply a sum of integrals of scalars, each 
of which is a contraction of the type appearing on the left side in (5). Using then (5) 
and (21) together in (28), and recalling that ~ on doAi, yields 



M 



dx 



M 



diM 



(29) 
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This gives rise to a covariant weak formulation of the problem: 

Find M e u + s.t. (i^(w),t;) = 0, \/ v e B2, (30) 

for suitable Banach spaces of functions Bi and B2, where the nonlinear weak form 
(F(-), •) can be written as: 



{F{u),v)^ g,j{A'^v^.^ + B'v') dx+ g.jC'v^ ds. (31) 

JM ' JdiM 

The notation {w, v) will represent the duality pairing of a function w in a Banach space 
B with a bounded linear functional (or form) w in the dual space B* . Depending on the 
particular function spaces involved, the pairing may be thought of as coinciding with 
the L^-inner-product through the Riesz Representation Theorem [ ] . The affine shift 
tensor u in (30) represents the essential or Dirichlet part of the boundary condition 
if there is one; the existence of u such that E = u\dg_M in the sense of the Trace 
operator is guaranteed by the Trace Theorem for Sobolev spaces on manifolds with 
boundary [4()], as long as E^ in (22) and doAi are smooth enough (see §2.4 below). 

The (Gateaux) linearization {DF(u)w,v) of the nonlinear form {F(u),v), 
necessary for both local solvability analysis and Newton-like numerical methods 
(of. [26]), is defined formally as: 



{DF{u)w,v) = ^{F{u + ew),v) 



(32) 



€=0 

This form is easily computed from most nonlinear forms {F(u),v) which arise from 
second order nonlinear elliptic problems, although the calculation can be tedious in 
some cases (as we will see shortly in the case of the constraints). 

The Banach spaces which arise naturally as solution spaces for the class of 
nonlinear elliptic systems in (30) are product spaces of the Sobolev spaces Wq'J^{M), 
or the related Besov spaces B^'P{A4). This is due to the fact that under suitable 
growth conditions on the nonlinearities in F, it can be shown (essentially by applying 
the Holder inequality) that there exists Pk,qk,rk satisfying 1 < Pk,q_k,Tk < 00 such 
that the choice 

Bi = iM)x---x W^p (M), B2 = Wl$ (Al) X . . . X {M), 

— + — = 1, rk>Tim\{pk,qk}, k = l,...,ne, (33) 
Pk qk 

ensures {F(u),v) in (31) remains finite for all arguments [47, 26]. 



2.4- The Sobolev Spaces W''-p{M) and H^{M) on Manifolds with Boundary 

For a type (r, s)-tensor T^j = 
satisfying |/| = r, jjj — s, define 



For a type (r, s)-tensor T^j — 7"°^"^ °bjb2 - b ' "^-^i^^re / and J are (tensor) multi-indices 



\T'j\^{T'jT\niLr""f'. (34) 
Here, 7/j and 7^'' are generated from the Riemannian d-metric ''jab on A4 as: 

lij = lablcd • • • 7p,, 1" = r'l'" ■■■Y'. (35) 

where |/| = | J| = m, producing m terms in each product. This is just an extension 
of the Euclidean P-norm for vectors in W^. For example, in the case of a 3-manifold, 
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taking |/| = 1, |J| = 0, 7a6 = Sab, gives \K\j\ = \K''\ = {K'^K^-tabf'^ = 

Employing the measure dx on M defined in (4), the L^-norm of a tensor on M., 
p £ [l,oo), is defined as: 

WT'jWlhm) =( [ I^Sr dx) , |lr^,|U»(A,) = ess sup \T'j{x)\. (36) 

The resulting i^'-spaces for 1 < p < oo are denoted as: 

LP{M) ^{T\j : llrSllLp(Ai) <^}- (37) 

Covariant (distributional) diS'erentiation of order m = \L\ (for some tensor multi- 
index L) using a connection generated by jab, or generated by possibly a different 
metric, will be denoted as either of: 

D"^T\j^T'j.j^, (38) 

where m should not be confused with a tensor index. The Sobolev semi-norm of a 
tensor is defined through (37) as: 

|L|=m 

and the Sobolev norm is subsequently defined as: 

\\T'j\\w-..iM)^ i E iT'jfw^^.iM)] ■ (40) 

yO<m<fe J 

The resulting Sobolev spaces of tensors are then defined as: 

W'''P{M) = {T\j : ||T^,||v^.,.(^) <oo }, 

Wl;'P{M) = { Completion of Cf^{M) w.r.t. || • Ww-.^iM) } , (41) 
where C^{Ai) is the space of C°°-tensors with compact support in A4. The space 
VFq '^(A^) is a special case of Wq'1){M), which can be characterized as: 

Wo${M) = {T\je W^^P : tr T^j.^ = on OqM, \L\ < k ^ I ] . (42) 

The spaces M^'^'P(A^) and W^o'^(A^) are separable (1 < p < oo) and refiexive 
{1 < p < oo) Banach spaces. The dual space of bounded linear functionals on 
W^''P{Ai) can be shown (in the sense of distributions, cf. [ ]) to be iy-'=^'?(7W), 
l/p+\/q = 1, which is itself a Banach space when equipped with the dual norm: 

\\f\\w-^HM)= sup „ ll^""^^ , - + - = l. (43) 

The Hilbert space special case of p = 2 is given a simplified notation: 

H''{M) = W''-^{M), H-''{M) = W-'''^{M), (44) 

with the same convention used for the various subspaces of H'^{A4) such as Hq{A4) 
and Hq The norm on H^{M.) defined above is induced by an L^-based inner- 

product as follows: \\T[j\\Hk(^M) = {T[j,T[jfll^j^y where 

{T'j,S[j)l2^m)= f T'jS\niLY^' dx, (45) 

JM 
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and where 

(T'j.S'j)h^(M)^ E (46) 

0<r?j.<fe 

The Banach spaces W^^'^ and their various subspaces satisfy various relations among 
themselves, with the L^-spaces, and with classical function spaces such as and C*^'". 
These relationships are characterized by a collection of results known as embedding, 
compactness, density, trace, and related theorems. A number of these are fundamental 
to approximation theory for elliptic equations, and will be recalled when needed below. 



2.5. Weak Formulation Example 

Before we derive a weak formulation of the Einstein constraints, let us consider a 
simple example to illustrate the idea. Here we assume the 3-metric to be flat so that 
V is the ordinary gradient operator and • the usual inner product. Let M. represent 
the unit sphere centered at the origin (with a single chart inherited from the canonical 
Cartesian coordinate system in M"^), and let dM. denote the boundary, satisfying the 
boundary assumption in equation (1). Consider now the following semilinear equation 
on M: 

- V • {a{x)Vu{x)) + 6(x, u{x)) = in M, (47) 
n{x) ■ {a{x)\Iu{x)) + c{x,u{x)) =Q on diM, (48) 
u{x) = f{x) on doM, (49) 

where n{x) : dAi t-^ M'^ is the unit normal to dAi, and where 

a : M >-^R^''^, b:MxR>-^R, (50) 
c : diM X M M, / : doM M. (51) 

To produce a weak formulation, we first multiply by a test function v G Hq d{A4) 
(the subspace of H^{A4) which vanishes on the Dirichlet portion of the boundary 
doM), producing: 

(-V • (aVu) + &(x,w))w = 0. (52) 

M 

After applying the flat space version of the divergence theorem, this becomes: 

/ (aVu) • Vw dx — / v{aS/u)-nds+ / b{x,u)vdx — 0. (53) 
Jm JdM Jm 

The boundary integral is reformulated using the boundary conditions as follows: 



v{aVu) ■ n ds — — / c{x,u)v ds. (54) 

dM JdiM 

If the boundary function / is regular enough so that / e H^/'^{doM), then from 
the Trace Theorem [32], there exists u G H^{Ai) such that / = uIq^m in the sense 
of the Trace operator. Employing such a function u G H^{A4), the weak formulation 
has the form: 

Find ueu + HI j^{M) s.t. {F{u),v) =0, V € HIo{M), (55) 
where from equations (53) and (54), the nonlinear form is defined as: 

{F{u),v) ^ / {aS/u ■ \/v + b{x,u)v) dx + c{x,u)v ds. (56) 

Jm JdiM 
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The "weak" formulation of the problem given by equation (55) imposes only one order 
of differentiability on the solution u, and only in the weak or distributional sense. 
Under suitable growth conditions on the nonlinearities h and c, it can be shown that 
this weak formulation makes sense, in that the form {F{-), ■) is finite for all arguments. 

To analyze linearization stability, or to apply a numerical algorithm such as 
Newton's method, we will need the bilinear linearization form {DF(u)w, v), produced 
as the formal Gateaux derivative of the nonlinear form {F(u),v): 



{DF{u)w,v) 



de 



{F{u + ew),v) 



e=0 



d_ 
de 



(aV(u + ew) • Vw + &(a:, M + ew)ti) da; + / c{x^u + ew)vds 
M JdiM 



'=0 



M 



a\7w ■ Vu - 



db{x, u) 



wv dx + 

ou J .IdiM 



dc{x, u) 
du 



wv ds. (57) 



Now that the nonlinear weak form {F(u),v) and the associated bilinear linearization 
form {DF(u)w,v) are defined as integrals, they can be evaluated using numerical 
quadrature to assemble a Galerkin-type discretization; this is described in some detail 
below in the case of a finite element-based Galerkin method. 



2. 6. Weak Formulation of the Constraints 

The Hamiltonian constraint (14) as well as the momentum constraint (15), taken 
separately or as a system, fall into the class of second-order divergence-form elliptic 
systems of tensor equations in (20)-(22). Therefore, we will follow the same plan as 
in §2.4 in order to produce the weak formulation (30)-(31). However, we now employ 
the conformal metric ^ab from the preceding section to define the volume element 
dx = y/dei ^ab dx^dx^dx^ and the corresponding boundary volume element ds, and 
for use as the manifold connection for covariant differentiation. The notation for 
covariant differentiation using the conformal connection will be denoted Da as in the 
previous section, and the various quantities from §2.4 will now be hatted to denote 
our use of the conformal metric. For example, the unit normal to dM. will now be 
denoted ff". 

Consider now the principle parts of the Hamiltonian and momentum constraint 
operators of the previous section: 

A(/. = L>a£'>, bb{LWY^ = bbib'^w'' + b'^w -\^''''b^w). (58) 

o 

Employing the covariant divergence theorem in equation (5) leads to covariant versions 
of the Green identities 

/ ^k(l)dx+ I (bacfyib^-tl^) dx= I haipb^fjids (59) 

JM Jm JdM 

and 

/ VabbiLWr^ dx + [ {LWr'-bbVadx^ [ hbVaiLWr'' ds, (60) 
Jm Jm Jom 

for smooth functions in C°°{M). These identities extend to W^'P{M) using a standard 
approximation argument, since C^{M.) is dense in W^'P{M.) (cf. Theorem 2.9 in [ "']). 
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Due to the symmetries of {LW)°'^ and ^"■^ , the second integrand in (60) can be 
rewritten in a completely symmetric form. To do so, we first borrow the hnear strain 
or (symmetrized) deformation operator from elasticity: 

(£;y)«6 ^ 1 (i)V" + b^V^) . (61) 



We can then write the operator {LWY^ in terms of {EW)""^ as follows: 

{LWy^ ^2n{EWY^ + \T^bcW'', with ^ = 1 and A = -^. (62) 

This makes it clear that the momentum constraint operator Df,(LT4^)°^ is a covariant 
version of the linear elasticity operator for a homogeneous isotropic material (cf. [ ]), 
for a particular choice of Lame constants. In particular, in the flat space case where 
7a{, = Sab, the operator becomes the usual linear elasticity operator: 

bbiLWy'' -^{2txeab{W) + Xe,,{W)5ab)^^, eab{W) ^]^{Wa^k + Wb^a) , (63) 

where Wa,b = dWa/dx^ denotes regular partial differentiation. Employing the 
operator {EW)""^ leads to a symmetric expression in the Green identity (60): 

{LWY'^bbVa = \ ({LWY'^bbVa + {LWT''baVb) (64) 

= {LWr\EV)ab (65) 

= 2^i{EWT\EV)ab + ^xr'b.w (^bbVa + baVb) (66) 

= 2fi{EWr\EV)ab + XbaWbbVK (67) 

While it is clear by inspection that the first operator in (58) is formally self-adjoint 
with respect to the covariant L^-inner-product defined in (45) , reversing the procedure 
in (67) implies that the same is true for the second operator in (58). In other words, 
the following holds (ignoring the boundary terms) : 

(A0, = (0, AV')l2(^), V(^,V', (68) 

{bbiLWr", F")l2(>i) = {W, bb{LVr')mM), V\ (69) 

To make it possible to write the Hamiltonian constraint and various related 
equations in a concise way, we now introduce the following nonlinear function 
P{<t>) = P((/), VF", a;''), where the explicit dependence on (and sometimes also the 
dependence on W"') is suppressed to simplify the notation: 

P{4>) = ^R4>^ + ^{trKf4>' + ^{*Aab + (LW^)„fc)20-6 + (70) 
The first and second functional derivatives with respect to are then as follows: 

P'{4>) =\R4>+^{irKf(l^^~\{*Aab + {LW)abfr' --i^pr^ (71) 

P"{4>) ^\r+ ^(trX)204 ^ 7^*4^^ ^ {LW)ab?4>-' + Q^pr^- (72) 

When working with the Hamiltonian constraint separately, we will often use these 
expressions involving P(-); when working with the coupled system, we will usually 
write the polynomial explicitly in order to indicate the coupling terms. 

We now take the inner product of (14) with a test function ip, which we assume to 
vanish on doM. (Again, d^M may have zero measure, or it may be all or only a piece 
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of dAi.) After use of the Green identity (59) and the Robin boundary condition (19) 
we obtain the following form (F^i^ct)) , tjj) (nonhnear in 0, but hnear in -ip) for use in 
the weak formulation in (30): 

(Fh((?!)),V') = / {c(l) - z)^; ds + I P'{(j))^p dx+ I bacfib"^ dx. (73) 

For the momentum constraint we take the inner product of (15) with respect to a 
test vector field V"" (again assumed to vanish on d^M.) and similarly use the Green 
identity (60) and the Robin condition (18) to obtain a form (Fm(V1^°), V'^) (linear in 
both W"" and V"" in this case) having the expression 



(Fm(W^'^),0 



diM 



{C\W^ -^Z'')Vads- 



M 



^"D^trK + Sirf K dx 



M 



(2pi{EWr''{EV)ab + XbaWb.V') dx, 



(74) 



where we have used (67). We will take /i = 1 and A = —2/3 in (74), but for the 
moment we will leave them unspecified for purposes of the discussion below. 

Ordering the Hamiltonian constraint first in the system (20), and defining the 
product metric Qij and the vectors and appearing in (27) and (31) as: 



1 

7ab 



(75) 



produces a single nonlinear weak form for the coupled constraints in the form required 
in (30), where 

{F{u),v) = (F([0, W-]), [V-, V-]) - (Fh(0), ^) + {Fm{W-), 
{[c(j) -z]ip+ [C^W'' - Z"] Va) ds+ [ (l(j)^b''tTK + SirjA Va dx 

^ dx 



M 



^-R4> + ^(tri^)205 - l(*i,, + {LW)abfr' - 2vr/5(^ 



-3 



M 



baCt^b^^J + 2fi{EWy''{EV)ab + XbaWbbV'') dx 



(76) 



While we have completely specified the weak form of the separate and coupled 
constraints on a manifold with boundary in a formal sense, they can be shown to 
be well-defined (and individually well-posed) in a more precise mathematical sense; 
see [10, 17, IN] for an analysis and a survey of the collection of existence, uniqueness, 
and stability results. For a particular situation, we must specify the particular 
combination of the boundary conditions (20)-(21) on a splitting of the boundary 
{dAi) into Dirichlet {doA4) and Robin {diA4) parts. This is quite problem dependent; 
in numerical simulation one typically computes solutions to the constraints in the 
interior of a large box or sphere. On the surface of the sphere one employs Robin 
and vector Robin conditions similar to those given in ['m], which fit the framework 
in (18) and (19). In addition, one often constructs black holes topologically by 
requiring the conformal metric to obey an isometry through one or more smaller 
non-overlapping spheres internal to the domain boundary. The isometry generates a 
boundary condition on the conformal factor which is well understood only when jab 
is flat. Even in this case, the exact corresponding boundary condition on is not 
known, but is likely to appear in the form of (18). Solvability of both constraints rests 
delicately on the boundary condition choices made. 
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2. 7. Gateaux Linearization of the Weak Formulation 



We will take the formal Gateaux-derivative of the nonlinear form in 
equation (76) above, to produce a linearization form for use in local solvability analysis 
through the Implicit Function Theorem, and for use in Newton-like iterative solution 
methods (cf. [ ]). Defining an arbitrary variation direction w = we compute 

the Gateaux-derivative of the nonlinear form as follows: 



d 
d~e 



diM 



([c(0 + ei)~z\^+ [C\{W'' + eX") - Z'^] K) ds 



e=0 



d 

d~e 



M 



d 

d~e 



+\Da{W'' + eX'')DbV^\ dx 

J c=0 

(f) + eS^fWtvK + Svrj" ) K dx 



M 



e=0 



d 

de 



M 



-^{*Aab + {L{W + eX))ab)^{4> + eCr' - 2ti p{cf> + ee)"^ j V dx 



(77) 



After some simple manipulations using the product and chain rules, we are left with 
the following bilinear form (for fixed [0, W^"]), linear separately in each of the variables 
[S^X"] and [V',^"]: 



diM 



{c(,i; + C\xWa) ds 



M 



DaiD"-ij + 2^i{EXriEV)ab + XDaX'^DbV' dx 



dx 



M 



{LXY''iP dx+ [ U(j)'^WtrK) Va^ dx. 



(78) 



Note that the first two volume integrals and the surface integral are completely 
symmetric in their arguments, and represent the symmetric part of the bilinear form. 
The third and fourth volume integrals are nonsymmetric in their arguments; the 
third volume integral represents the linearized coupling of W"' into the Hamiltonian 
constraint, and the fourth volume integral represents the linearized coupling of the 
conformal factor (h into the momentum constraint. 
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2.8. Weak Formulations Arising from Energy Functionals 

Due to the fact that the principle parts of the Hamiltonian and momentum operators 
produced by the conformal decomposition are self-adjoint, the weak formulations 
individually arise naturally as the Euler conditions for stationarity of associated 
(energy) functionals. It is straight-forward to verify that the following energy 
functionals 

M^) = \ I c(0 - z)c^ ds+ f P{^, W) dx+l f bacj^b'^c^ dx, (79) 



2 JdiM JM 2 jj^ 



JuiW) = \ [ {C%W' - Z'^) Wads+ I (l4,''WtYK + SnjA Wa dx (80) 

^ JdiM JM V'J / 



M 



^(2fi{EWr\EW)ab + XbaWbbW''^ dx, (81) 

each separately give rise to the individual weak forms (73) and (74), respectively. One 
computes the Gateaux derivative of Jh(0, W^") with respect to 0, and the Gateaux 
derivative of Jm(0, W"") with respect to W"^, and then sets them to zero: 

^JH((/. + eV) =(J{j(0),V^) = (Fh(0),V) = O, (82) 



de 



= {MW). V'') = {Fm{W'),V'') = 0. (83) 

This discussion is only formal, but can be made rigorous. 

Unfortunately, while the individual constraints each arise as Euler conditions for 
stationarity of the separate energy functionals above, the coupled constraints do not 
arise in this way from any coupled energy. This follows easily from the fact that the 
combined linearization bilinear form in (78) is not symmetric. This can also be verified 
directly by considering the most general possible expression for the total energy: 

JtotaM, = M4>, W) + Jm(0, W) + Jr(0, W), (84) 

where Jy{{4>,W°-) and Jm(<^, VF") are as defined above, and where Jr,((/), ly^) is 
the remainder term in the energy which must account for the coupling terms in 
the combined weak form (76). It is easy to verify that the Euler condition for 
stationarity places separate conditions on the Gateaux derivative of Jr(0, W"') which 
are impossible to meet simultaneously. 

This lack of a variational principle for the coupled constraints limits the number of 
techniques available for analyzing solvability of the coupled system; the existing near- 
CMC results [7, 16] and the non-CMC (far-from-CMC) results [17, IS] are actually 
based on fixed-point arguments; however, variational arguments are used to solve 
the individual constaint equations in [16, 17, 18] as part of the overall fixed-point 
argument. 



3. Adaptive Finite Element Methods (AFEM) 

In this section we give a brief description of Galerkin methods, finite element methods, 
and adaptive techniques for covariant nonlinear elliptic systems. We also derive a 
posteriori error indicators for driving adaptivity, and finish the section with some a 
priori error estimates for general Galerkin approximations to the Hamiltonian and 
momentum constraints. Expanded versions of this material, including proofs of all 
results, can be found in [26, 27, 28]. 
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3.1. Petrov-Galerkin Methods, Galerkin Methods, and Finite Element Methods 

A Petrov-Galerkin approximation of the solution to (30) is the solution to the following 
subspace problem: 

Find {uh - Uh) e Uh C Bi s.t. {F{uh),v) = 0, (85) 

for some chosen subspaces Uh and Vh, where (imi{Uh) — dim(V/i) — n, and where 
the discrete Dirichlet function Uh approximates u (e.g. an interpolant). A Galerkin 
approximation refers to the case that Uh — Vh. 

A finite element method is simply a Petrov-Galerkin or Galerkin method in which 
the subspaces Uh and Vh are chosen to have the extremely simple form of continuous 
piecewise polynomials with local support, defined over a disjoint covering of the domain 
manifold M by elements. For example, in the case of continuous piecewise linear 
polynomials defined over a disjoint covering with 2- or 3-simplices (cf. Figure 3.1), 
the basis functions are easily defined element-wise using the unit 2-simplex (triangle) 
and unit 3-simplex (tetrahedron) as follows: 



<Po{x,y) = 1- 
<^i(^,y) = i 

4>2{x,y) = y 



(f^o{S:,y,z) = 1-x-y 

(pi{x,y,z) = y 

(h{x,y,z) = X 

4>3{x,y,z) = z 



Global basis functions are then defined, as in the right-most picture in Figure 3.1, 



Figure 1. Reference and arbitrary 2- and 3-simplex elements, and a global {2D) 
basis function. 



by combining the support regions around a given vertex, and extending the unit 
simplex basis functions to each arbitrary simplex using coordinate transformations. 
If the manifold domain can be triangulated exactly with simplex elements, then the 
coordinate transformations are simply affine transformations. Note that in this sense, 
finite element methods are by their very nature defined in a chart-wise manner. 
Quadratic and high-order basis functions are defined analogously. 

The above basis functions clearly do not form any subspace of C^(A^), the space 
of twice continuously differentiable functions on M, which is the natural function 
space in which to look for the solutions to second order elliptic equations. This is due 
to the fact that they are discontinuous along simplex faces and simplex vertices in the 
disjoint simplex covering of M. However, one can show [ ] that in fact: 

T4 = span{0i,...,0„}c<5(X), MdR^, (86) 
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so that these continuous, piecewise defined, low-order polynomial spaces do in fact 
form a subspace of the solution space to the weak formulation of the class of second 
order elliptic equations of interest. Making then the choice Uh = span{(/>i, (f)2^ ■ ■ ■ , 0n}, 
Vh = span{-0i, -02, . . . , '0n}i equation (85) reduces to a set of n nonlinear algebraic 
relations (implicitly defined) for the n coefficients {aj} in the expansion 

n 

Uh ^ Uh + ^aj<j>j. (87) 

In particular, regardless of the complexity of the form {F{u),v), as long as we can 
evaluate it for given arguments u and v, then we can evaluate the nonlinear discrete 
residual of the finite element approximation Uh as: 

n 

{F{uh + ^aj(j)j),ijji), z = l,...,n. (88) 

Since the form {F{u), v) involves an integral in this setting, if we employ quadrature 
then we can simply sample the integrand at quadrature points; this is a standard 
technique in finite element technology. Given the local support nature of the functions 
(jjj and Tpi, all but a small constant number of terms in the sum ^re zero 

at a particular spatial point in the domain, so that the residual is inexpensive to 
evaluate when quadrature is employed. 

The two primary issues in applying this approximation method are then: 

(i) The approximation error ||u — it^Hx, for various norms X, and 

(ii) The computational complexity of solving the n algebraic equations. 

The first of these issues represents the core of finite element approximation theory, 
which itself rests on the results of classical approximation theory. Classical references 
to both topics include [51, 52, 53]. The second issue is addressed by the complexity 
theory of direct and iterative solution methods for sparse systems of linear and 
nonlinear algebraic equations, cf. [54, 55], and by the use of adaptive techniques to 
minimize the size n of the discrete space that must be constructed to reach a specific 
approximation quality. 

3.2. A Priori Error Estimates for the Gonstraint Equations 

We first outline an approximation result for general Galerkin approximations to 
solutions of the momentum constraint. It is referred to as a quasi-optimal error 
estimate, in that it establishes a bound on the error jju — that is within a 

constant of being the error in the best possible approximation. 

To understand this result, we begin with the two (Hilbert vector) spaces H and 
V, where in our setting of the momentum constraint, we have H = L^{A4) and 
V = Hq We will stay with the abstract notation involving H and V for clarity. 

The weak form of the momentum constraint can be shown to have the following form: 

Find u&V s.t. A{u, v) = F{v), Vw £ V, (89) 
where the bilinear form A(u, v) : V x V t-^ Wis bounded 

A{u,v) < M\\u\\v\\v\\v, yu,veV, (90) 
and V-coercive (satisfying a Carding inequality): 

m\\u\\^ < K\\u\\% + A{u,u), Wu&V, where m > 0, (91) 
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and where the hnear functional F{v) : V ^ is bounded and thus hes in the dual 
space V*: 

F{v) < L\\v\\v, yv e V. 

It can be shown that the weak formulation of the momentum constraint fits into this 
framework; to simplify the discussion, we have assumed that any Dirichlet function u 
has been absorbed into the linear functional F(v) in the obvious way. Our discussion 
can be easily modified to include approximation of u hy Uh- 

Now, we are interested in the quality of a Galerkin approximation: 

Find Uh e Vh C V s.t. A{uh,v)^A{u,v)^F{v), \fv £ Vh C V. (92) 

We will assume that there exists a sequence of approximation subspaces Vh C V 
parameterized by h, with Vh-^ C Vh^ when /i2 < ft-i, and that there exists a sequence 
{a/i}, with lim/i^o o.h — 0, such that 

\\u — Uh\\H < — u/illi/, when A(u — u^, v) — 0, Vu G V/i C V. (93) 

The assumption (93) is very natural; in our setting, it is the assumption that the 
error in the approximation converges to zero more quickly in the L^-norm than in the 
H^-norm. This is easily verified in the setting of piecewise polynomial approximation 
spaces, under very mild smoothness requirements on the solution u. Under these 
assumptions, we have the following a priori error estimate. 

Theorem 3.1 For h sujjiciently small, there exists a unique approximation Uh 
satisfying (92), for which the following quasi-optimal a priori error hounds hold: 

\\u-Uh\\v <C mi \\u-v\\v, (94) 
veVh 

\\u-Uh\\H <Cah inf \\u-v\\v, (95) 

where C is a constant independent of h. If K < in (91), then the above holds for all 
h. 

Proof See [26, 27, 28]; also [56]. □ 

As we did previously for the momentum constraint, we now outline an 
approximation result for general Galerkin approximations to solutions of the 
Hamiltonian constraint. Again, it is referred to as a quasi-optimal error estimate, 
in that it establishes a bound on the error \\u — Uh\\x that is within a constant of 
being the error in the best possible approximation. 

We begin again with the two Hilbert spaces H and V, where again we have 
H = L'^{M) and V = II^jj{A4). We are given the following nonlinear variational 
problem: 

Find u eV s.t. A{u,v) + {B{u),v) ^ F{v), \fv e V, (96) 
where the bilinear form A(u, v) : V x V t-^ R is bounded 

A{u,v) < M\\u\\v\\v\\v, yu,veV, (97) 

and V-eUiptic: 

'Tijlujly < A(u,w), '^ueV, where to > 0, (98) 

where the linear functional F(v) : V R is bounded and thus lies in the dual space 
V*: 

F{v) < L\\v\\v, yv e V, 
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and where the nonhnear form {B{u),v) : V x V t-^ M. is assumed to be monotonic: 

< {B{u) - B{v),u-v), yu,veV, (99) 

where we have used the notation: 

{B{u) - B{v),w) = {B{u),w) - {B{v),w). (100) 

We are interested in the quahty of a Galerkin approximation: 

Find Uh e Vh s.t. A{uh,v) + {B{uh),v)^F{v), VveVh, (101) 

where Vh C V. We will assume that {B{u),v) is bounded in the following weak sense: 
li u G V satisfies (96), if Uh S Vh satisfies (101), and if w e Vh, then there exists a 
constant K > such that: 

{B{u)- B{uh),u-v) <K\\u-Uh\\v\\u-v\\v. (102) 

It can be shown that the weak formulation of the Hamiltonian constraint fits into this 
framework. We have again assumed that any Dirichlet function u has been absorbed 
into the various forms in the obvious way to simplify the discussion. The discussion 
can be modified to include approximation of u by Uh- 

Again, we are interested in the quality of a Galerkin approximation Uh 
satisfying (101), or equivalently: 

A{u - Uh, v) + {B{u) - B{uh),v) ^0,yveVhC V. 

As before, we will assume that there exists a sequence of approximation subspaces 
Vh <Z V parameterized by h, with Vh-^ C Vh2 when /i2 < hi, and that there exists a 
sequence {ah}, with lim/i^o i/i = 0, such that 

\\u-Uh\\H <ah\\u-Uh\\v, (103) 

holds whenever Uh satisfies (101). The assumption (103) is again very natural; see 
the discussion above following (93). Under these assumptions, we have the following 
a priori error estimate. 

Theorem 3.2 The approximation Uh satisfying (101) obeys the following quasi- 
optimal a priori error bounds: 

\\u-uh\\v < C inf \\u-v\\v, (104) 
\\u - Uh\\H < Cah inf (105) 

where C is a constant independent of h. 

Proof See [26, 27, 28]. □ 

3.3. Adaptive Finite Element Methods (AFEM) 

A priori error analysis for the finite element method for addressing the first issue is 
now a well-understood subject [51, 57]. Much activity has recently been centered 
around a posteriori error estimation and its use in adaptive mesh refinement 
algorithms [ nS, 59, 60, 61, 62, 63]. These estimators include weak and strong residual- 
based estimators [59, 60, 61], as well as estimators based on the solution of local 
problems [64, 65]. The challenge for a numerical method is to be as efficient as 
possible, and a posteriori estimates are a basic tool in deciding which parts of the 
solution require additional attention. While the majority of the work on a posteriori 
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estimates has been for linear problems, nonlinear extensions are possible through 
linearization theorems (cf. [61, G2] and the discussion of the error estimator employed 
by FETK later in this paper). The solve-estimate-refine structure in simplex-based 
adaptive finite element codes such as FETK [2()] and PLTMG [oS], exploiting these 
a posteriori estimators, is as follows: 
Algorithm 1 (Adaptive multilevel finite elements) 

• While {\\u — Uh\\x > e) do: 

(i) Find (uh — Uh) £ Uh C Bi such that: 

{F{uh),v)=o, \/veVhCB2. 

(ii) Estimate \\u — Uh\\x over each element. 

(Hi) Initialize two temporary simplex lists as empty: Ql — Q2 = 0. 

(iv) Place simplices with large error on the "refinement" list Ql. 

(v) Bisect all simplices in Ql (removing from Ql), and place any nonconforming 
simplices created on the list Q2. 

(vi) Ql is now empty; set Ql = Q2, Q2 = 0. 

(vii) If Ql is not empty, goto (5). 

• End While. 

The conformity loop (5)-(7), required to produce a globally "conforming" mesh 
(described below) at the end of a refinement step, is guaranteed to terminate in a 
finite number of steps (cf. [66, 67]), so that the refinements remain local. Element 
shape is crucial for approximation quality; the bisection procedure in step (5) is 
guaranteed to produce nondegenerate families if the longest edge is bisected in 
two dimensions [08, 69], and if marking or homogeneity methods are used in three 
dimensions [22, 23, 70, 71, 72, 73]. Whether longest edge bisection is nondegenerate 
in three dimensions apparently remains an open question. Figure 2 shows a single 
subdivision of a 2-simplex or a 3-simplex using either 4-section (left-most figure), 8- 
section (fourth figure from the left), or bisection (third figure from the left, and the 
right-most figure). The paired triangle in the 2-simplcx case of Figure 2 illustrates 



Subdividing 2-slmpllces (tdangles) 



Subdividing 3-slmpllces (tetrahedra) 




Figure 2. 

bisection. 



Refinement of 2- and 3-simpIices using 4-section, 8-section, and 



the nature of conformity and its violation during refinement. A globally conforming 
simplex mesh is defined as a collection of simplices which meet only at vertices and 
faces; for example, removing the dotted bisection in the third group from the left in 
Figure 2 produces a non-conforming mesh. Non-conforming simplex meshes create 
several theoretical as well as practical implementation difficulties, and the algorithms 
in FETK (as well as those in PLTMG [')^] and similar simplex-based adaptive 
codes [26, 23, 74, 75, 76]) enforce conformity using the above queue swapping strategy 
or a similar approach. 

Addressing the complexity of Step 1 of the algorithm above, Newton methods are 
often the most effective: 
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Algorithm 2 (Damped-inexact-Newton) 

• Let an initial approximation u be given. 

• While (\{F{u),v)\ > e for any v) do: 

(i) Find w such that: 

{DF{u)w, v) = -{F(u), v} + r, V v. 
(ii) Set u = u + \w. 

• End While. 

The bilinear form {DF(u)w,v) in the algorithm above is simply the (Gateaux) 
linearization of the nonlinear form {F(u),v), defined formally as: 



{DF{u)w,v) = ^{F{u + ew),v) 



(106) 



This form is easily computed from most nonlinear forms {F(u),v) which arise from 
second order nonlinear elliptic problems, although the calculation can be tedious in 
some cases (as in the case of the constraints in general relativity) . The possibly nonzero 
"residual" term r is to allow for inexactness in the Jacobian solve for efficiency, which 
is quite effective in many cases (cf. [77, 7H, 79]). The parameter A brings robustness 
to the algorithm [79, 80, SI]. If folds or bifurcations are present, then the iteration is 
modified to incorporate path- following [82, 83]. 

As was the case for the nonlinear residual {F{-), •), the matrix representing the 
bilinear form in the Newton iteration is easily assembled, regardless of the complexity 
of the bilinear form {DF{-)-, •). In particular, the algebraic system for w = X]J=i Pj'Pj 
has the form: 



AU = F, U,^f3„ (107) 



where 



A,, = {DF{uh + ^afc(/)fe)0j, (108) 
fc=i 

n 

F, ={F{uh + Y,aj(j),),^,). (109) 
i=i 

As long as the integral forms {F{-), •) and {DF{-)-, •) can be evaluated at individual 
points in the domain, then quadrature can be used to build the Newton equations, 
regardless of the complexity of the forms. This is one of the most powerful features 
of the finite element method, and is exploited to an extreme in the code FETK (see 
Section 6 and [2()]). It should be noted that there is a subtle difference between 
the approach outlined here (typical for a nonlinear finite element approximation) and 
that usually taken when applying a Newton-iteration to a nonlinear finite difference 
approximation. In particular, in the finite difference setting the discrete equations are 
linearized explicitly by computing the Jacobian of the system of nonlinear algebraic 
equations. In the finite element setting, the commutativity of linearization and 
discretization is exploited; the Newton iteration is actually performed in function 
space, with discretization occurring "at the last moment" in Algorithm 2 above. 

It can be shown that the Newton iteration above is dominated by the 
computational complexity of solving the n linear algebraic equations in each iteration 
(cf. [77, 84]). Multilevel methods are the only known provably optimal or nearly 
optimal methods for solving these types of linear algebraic equations resulting 
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from discretizations of a large class of general linear elliptic problems [84, 85, 86]. 
An obstacle to applying multilevel methods to the constraint equations in general 
relativity and to similar equations is the presence of geometrically complex domains. 
The need to accurately represent complicated domain features and boundaries with an 
adapted mesh requires the use of very fine mesh simply to describe the complexities of 
the domain. This may preclude the use of the solve-estimate-refine structure outlined 
above in some cases, which requires starting with a coarse mesh in order to build 
the approximation and linear algebra hierarchies as the problem is solved adaptively. 
In this situation, algebraic or agglomeration/aggregation-based multilevel methods 
can be employed [87, 88, 89, 90, 91, 92, 93, 94, 95]. A fully unstructured algebraic 
multilevel approach is taken in FETK, even when the refinement hierarchy is present; 
see Section 6 below and also [ ] for a more detailed description. 

5.^. Residual-Based A Posteriori Error Indicators 

There are several approaches to adaptive error control, although the approaches based 
on a posteriori error estimation are usually the most effective and most general. While 
most existing work on a posteriori estimates has been for linear problems, extensions 
to the nonlinear case can be made through linearization. To describe one such result 
from [li()], we assume that the d-manifold M has been exactly triangulated with a 
set iS of shape-regular d-simplices (the finite dimension d is arbitrary throughout this 
discussion). A family of simplices will be referred to here as shape-regular if for all 
simplices in the family the ratio of the diameter of the circumscribing sphere to that 
of the inscribing sphere is bounded by an absolute fixed constant, independent of the 
numbers and sizes of the simplices that may be generated through refinements. (For 
a more careful definition of shape-regularity and related concepts, see [51].) It will be 
convenient to introduce the following notation: 



S = Set of shape-regular simplices triangulating M 

N{s) = Union of faces in simplex set s lying on On-M 

T{s) = Union of faces in simplex set s not in A/'(s) 

I^{s) = A/'(.s)UX(s) 

= [j {seS \ sCls^^, where s eS } 

LOf = U { Se5 I /nS^0, where/e.F} 

hs = Diameter (inscribing sphere) of the simplex s 

hf = Diameter (inscribing sphere) of the face /. 



When the argument to one of the face set functions Af, T, or !F is in fact the entire 
set of simplices S, we will leave off the explicit dependence on S without danger of 
confusion. Referring forward briefiy to Figure 5 will be convenient. The two darkened 
triangles in the left picture in Figure 5 represents the set Wf for the face / shared by 
the two triangles. The clear triangles in the right picture in Figure 5 represents the 
set Ws for the darkened triangle s in the center (the set Ws also includes the darkened 
triangle). 

Finally, we will also need some notation to represent discontinuous jumps in 
function values across faces interior to the triangulation. To begin, for any face / e J\f, 
let Uf denote the unit outward normal; for any face / € X, take rif to be an arbitrary 
(but fixed) choice of one of the two possible face normal orientations. Now, for any 
V G L^{A4) such that v e C^{s) Vs G S, define the jump function: 

[v]f{x) ^ lim v{x + enf)— lim v{x — enf). 

e— >0+ e— tO~ 
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By analyzing the clement- wise volume and surface integrals in (31) and using some 
technical results on interpolation of functions, the following fairly standard result is 
derived in [ ' ]: 

Theorem 3.3 Let u e iy^'''(A^) he a regular solution of (20)-(22), or equivalently 
of (30)-(31), where some additional minimal assumptions hold as described 
in [26]. Then the following a posteriori error estimate holds for a Petrov-Galerkin 
approximation uu satisfying (85): 

\\u-Uh\\w^..^M)<c{y^Ti^\ , (110) 



where 



C = 2.max{C„C/} - maxli^y^i?)/'} • \\DF{u)-^\\ciw-^^^.w'^^), 



and where the element-wise error indicator rjg is defined as 



/.?||i?^-A»J|^,(,) + ^ hs\\[A^^n,]f\\l,^f) (111) 



2 



1/p 

Proof. See [26, 27, 28]. □ 

3. 5. An A Posteriori Error Indicator for the Constraints 

Here, we using the general indicator above, we can instantiate an estimator specifically 
for the constraints in general relativity. The Ph.D. thesis of Mukherjee [2.H] contains 
a residual-based error estimator for the Hamiltonian constraint that is equivalent to 
our estimator when the momentum constraint is not involved, in the specific case of 
p = q = r = 2. We consider first the Hamiltonian constraint, which can be thought of 
as an equation of the form (20)-(22). The error indicator from above now takes the 
form: 



fei{s} /eAA(s) / 

The momentum constraint also has the form (20)-(22), and the error indicator in [26] 
takes the form: 

vf=(h':\\l<l>'D''trK + 8nf-D,{LW,rX^,^^^ + ^- ^ hf\\ [mLW^T"] Wl,^^^ 

+ hf\\c^,wi:~z^ + h,{Lw,rXHf)] ■ (113) 
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Finally, the error indicator we employ for the coupled system is the P-weighted average 
of the two estimators above: 

^ {wH{vfY + WMivfyf' , (114) 

where the weights wh and wm satisfying wh + wm = 1 are determined heuristically. 

Note that in this special case the weighted sum estimator in (114) can be derived 
from the general system estimator in [2G] by defining the product space metric 
described in [ ] as follows: 



Wh 
















WMQab 


y = 









(115) 



and employing the coupled system framework from [26, 27]. However, in general 
different components of a coupled system could lie in different function spaces, and the 
norm appearing in the estimator in [2(), 27] would need to be modified. Equivalently, 
an estimator built from a weighted sum of estimators for individual components could 
be used as above. 



3.6. Gonvergence and Optimal Gomplexity of AFEM for the Gonstraints 

The following convergence result for AFEM applied the a general class of nonlinear 
elliptic equations that includes the Hamiltonian constraint appears in [28]. More 
general results which apply to the coupled system appear in [',]. Below, the energy 
norm is used as the natural norm for the analysis: 

|||u||| — A{u, u)^, 

where the bilinear form A(u,v) is defined in Section 2.6. 

Theorem 3.4 (Contraction) Let {Pk, Sk,Uk}k>o be the sequence of finite element 
meshes, spaces, and solutions, respectively, produced by AFEM. Let the initial mesh 
size Hq be sufficiently small so that a quasi-orthogonality inequality holds (see [28]) 
holds for {Pq, Sq, C/o}/£>o- Then, there exist constants 7 > and a G (0, 1), depending 
only on some freely specifiable AFEM parameters and the shape-regularity of the initial 
triangulation Pq, such that 

\\\u - Uk+if + Wk+i < « (Ilk - c^felll' + Wk) ■ 

Proof See [28]. □ 

The strict contraction result above makes possible the following complexity result 
from [28], which guarantees optimal complexity of the AFEM iteration for the 
Hamiltonian constraint equation under reasonable assumptions on the problem data. 
The approximation class As is precisely characterized in [27, 28]. 

Theorem 3.5 (Optimality) If data lies in approximation class Ag, then there exists 
a constant C such that 

\\\u - Uk+if + oscfc < C i#Pk - #Po)"' • 



Proof See [28]. □ 
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4. Fast Solvers and Preconditioners for AFEM 

We first introduce a superscript (j) to indicate the level of the stiffness matrix (or 
sometimes, the discretization operator) in the algebraic system (107). At certain 
places, we will drop the superscript for simplicity. The FE solution is sampled on 
the nodes, hence the number of internal nodes is equal to the number of degrees of 
freedom (DOF) or unknowns in the system (107). The total number of DOF on the 
finest level J is denoted by Nj = N . 

The stiffness matrix A'^^^ is ill-conditioned, with the condition number k,{A) of 
the system (107) growing 0(2^^ ) as — > oo (in the case of second order elliptic PDE). 
It is imperative to improve the condition of (107) by transforming the system to an 
equivalent one, namely by preconditioning 

where k(C'-'^^^^ A^^^ C^-'-'^^^) <C k{A'^^^). Moreover, the preconditioning matrix C^-'' 
has to be positive definite, and in some sense simple. One possible way to implement 
the above strategy is to determine a positive definite matrix B'^^\ having the following 
two properties: 

can efficiently be computed (usually, 0{Nj) is the desirable bound for the 
number of arithmetical operations when solving a linear system with coefficient 
matrix B^^'>), 

• A'-^^ and i?'--'^ are "almost" spectrally equivalent, i.e. 

Xbx^b'^^^x < A^'^\ < Abx^b'-^^x, X e 

with two positive constants As , with a small ratio . 

Since K^B^^r^^^ A^^'> B^^^^^'^) < then C^^'> = Sf^)"^ will be a good 

preconditioner choice. 

Solution of the algebraic system (107) by iterative methods has been the subject of 
intensive research because of the enormous practical impact on a number of application 
areas in computational science. For quality approximation in physical simulation, 
one is required to use meshes containing very large numbers of simplices leading 
to approximation spaces Sj with very large dimension Nj. Only iterative methods 
which scale well with Nj can be used effectively, which usually leads to the use of 
multilevel-type iterative methods and preconditioners. Even with the use of such 
optimal methods for (107), which means methods which scale linearly with Nj in 
both memory and computational complexity, the approximation quality requirements 
on Sj often force Nj to be so large that only parallel computing techniques can be 
used to solve (107). 

To overcome this difficulty one employs adaptive methods, which involves the 
use of a posteriori error estimation to drive local mesh refinement algorithms. This 
approach leads to approximation spaces Sj which are adapted to the particular target 
function u of interest, and as a result can achieve a desired approximation quality with 
much smaller approximation space dimension Nj than non-adaptive methods. One 
still must solve the algebraic system (107), but unfortunately most of the available 
multilevel methods and preconditioners are no longer optimal, in either memory or 
computational complexity. This is due to the fact that in the local refinement setting. 
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the approximation spaces Sj do not increase in dimension geometrically as they do 
in the uniform refinement setting. As a result for example, a single multilevel V- 
cycle no longer has linear complexity, and the same difficulty is encountered by other 
multilevel methods. Moreover, storage of the discretization matrices and vectors for 
each approximation space, required for assembling V-cycle and similar iterations, no 
longer has linear memory complexity. 

A partial solution to the problem with multilevel methods in the local refinement 
setting is provided by the hierarchical basis (HB) method [!)G, 97, 98]. This 
method is based on a direct or hierarchical decomposition of the approximation 
spaces Sj rather than the overlapping decomposition employed by the multigrid and 
Bramble-Pasciak-Xu (BPX) [99] method, and therefore by construction has linear 
memory complexity as well as linear computational complexity for a single V-cycle- 
like iteration. Unfortunately, the HB condition number is not uniformly bounded, 
leading to worse than linear overall computational complexity. While the condition 
number growth is slow (logarithmic) in two dimensions, it is quite rapid (geometric) 
in three dimensions, making it ineffective in the 3D local refinement setting. Recent 
alternatives to the HB method, including both BPX-like methods [100, 99] and 
wavelet-like stabilizations of the HB methods [101], provide a final solution to the 
condition number growth problem. It was shown in [102] that the BPX preconditioner 
has uniformly bounded condition number for certain classes of locally refined meshes 
in two dimensions, and more recently in [103, 104, 105, 106] it was shown that the 
condition number remains uniformly bounded for certain classes of locally refined 
meshes in three spatial dimensions. In [103, 105, 107], it was also shown that wavelet- 
stabilizations of the HB method give rise to uniformly bounded conditions numbers 
for certain classes of local mesh refinement in both the two- and three-dimensional 
settings. 

Ji^.l. Preliminaries on Optimal Preconditioners 

In the uniform refinement setting, the parallelized or additive version of the multigrid 
method, also known as the BPX preconditioner, is defined as follows: 

Xz.:=5]2^('^-2)^(^^^p))^p)^ ^^5^, (llg) 

j=0 i=l 

Only in the presence of a geometric increase in the number of DOF, the same 
assumption for optimality of a single classical (i.e. smoother acting on all DOF) 
multigrid or BPX iteration, does the cost per iteration remain optimal. In the 
case of local refinement, the BPX preconditioner (116) (usually known as additive 
multigrid) can easily be suboptimal because of the suboptimal cost per iteration. If the 
smoother is restricted to the space generated by fine or newly created basis functions, 
i.e. Sj := {Ij — Ij-i) Sj, then (116) corresponds to the additive HB preconditioner 
in [98]: 

J 

X^^u^Y.^^^'-'^ E ueSj. (117) 

j=0 i=Nj-i + l 

However, the HB preconditioner suffers from a suboptimal iteration count. Namely, 
if the algebraic system (107) is preconditioned by the HB preconditioner, the arising 
condition number is 0{J^) and 0(2"') in 2D and 3D, respectively. 
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In the local refinement setting, in order to maintain optimal overall computational 
complexity, the remedy is to restrict the smoother to a local space Sj which is typically 
slightly larger than the one generated by basis functions corresponding to fine DOF: 

{1,-1,-,) SjQSjCSj, (118) 

where /_,■ : L2{^) — > Sj denotes the finite element interpolation operator. The 
subspace generated by the nodal basis functions corresponding to fine or newly 
created degrees of freedom (DOF) on level j is denoted by {Ij — Ij~i) Sj. Nodes- 
equivalently, DOF-corresponding to Sj and their cardinality will be denoted by A/} 
and Nj , respectively. The above deficiencies of the BPX and HB preconditioners can 
be overcome by restricting the (smoothing) operations to Sj. This leads us to define 
the BPX preconditioner for the local refinement setting as: 

Xu ■.^Y.V^^-^) J2 u e Sj. (119) 

j"=o igXAj 

Remark 4.1 In order to prove optimal results on convergence, the basic theoretical 
restriction on the refinement procedure is that the refinement regions from each level 
forms a nested sequence. Let U,j denote the refinement region, namely, the union of 
the supports of basis Junctions which are introduced at level j . Due to nested refinement 
flj C ^j-i. Then the following nested hierarchy holds: 

fij C flj-i C ■ ■ ■ C flo ^ ft. 

Simply, the restriction indicates that tetrahedra of level j which are not candidates for 
further refinement will never be touched in the future. In practice, depending on the 
situation, the above nestedness restriction may or may not be enforced. We enforce 
the restriction in Lemma 4.1. In realistic situations, it is typically not enforced. 

4.2. Matrix Representations and Local Smoothing 

We describe how to construct the matrix representation of the preconditioners under 
consideration. Let the prolongation operator from level j — 1 to j be denoted by 

and also denote the prolongation operator from level j to J as: 

P, = Pf = . . . e M^■^><^^ 

where Pj is defined to be the rectangular identity matrix / G M.^-'^^-'-'^ . Then the 
matrix representation of (116) becomes [1()(S]: 

X = ^2J"('*-2)p^.p^t^ 

3=0 

One can also introduce a version with an explicit smoother Gj : 

X = Y.P,G,Pl 
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Throughout this article, the smoother Gj G M ^ is a symmetric Gauss-Seidel 
iteration. Namely, Gj = {Dj + Uj)^^ Dj{Dj + Lj)^^ where Aj = Dj + Lj + Uj with 
No = Nq. 

The matrix representation of (117) is formed from matrices Hj which are simply 
the tails of the Pj corresponding to newly introduced DOF in the fine space. In 
other words, Hj £ R^J^^i^j-^j-i) is given by only keeping the fine columns (the last 
Nj — Nj_i columns of Pj). Hence, the matrix representation of (117) becomes: 

XHB = X^2^('^-^)i?,i?j. 

3=0 

If the sum over i in (116) is restricted only to those nodal basis functions with 
supports that intersect the refinement region [109, 100, 102, 110], then we obtain the 
set called as onering of fine DOF. Namely, the set which contains fine DOF and their 
immediate neighboring coarse DOF: 

ONERING^^'^ = {onering{ii) : ii = Nj^i + 1, . . . , Nj}, 

where cmering(ii) = {ii, father s{ii)}. Now, the generic preconditioner (119) for local 
refinement transforms into the following preconditioner: 

Xu = X^2^('^-2) (u,cjJf^)cl^\ ueSj. (120) 

There are three popular choices for Afj . We outline possible BPX choices by the 
following DOF corresponding to: 

• (DOF-1) The basis functions with supports that intersect the refinement region 
fij [109, 100, 102]. We call this set onering^-''). 

• (DOF-2) The basis functions with supports that are contained in [ ■. . ]■ 

• (DOF-3) Created by red refinement and their corresponding coarse DOF. 

Here, red refinement refers to quadrasection or octasection in 2D and 3D, respectively. 
Green refinement simply refers to bisection. 

The interesting ones are DOF-1 and DOF-3 and we would like to elaborate on 
these. In the numerical experiments reported in [i ], DOF-1 was used. For the 
provably optimal computational complexity result in Lemma 4.1 DOF-3 is used. 

4-2.1. The Sets DOF-1, DOF-3 and Local Smoothing Computational Complexity The 
set DOF-1 can be directly determined by the sparsity pattern of the fine-fine subblock 
A22 of the stiffness matrix in (123). Then, the set of DOF over which the BPX 
preconditioner (120) smooths is simply the union of the column locations of nonzero 
entries corresponding to fine DOF. Using this observation, HB smoother can easily be 
modified to be a BPX smoother. 

DOF-3 is equivalent to the following set: 

AG- - - iV,-i + 1, . . . , N,} : ^ z - 1, . . . , iV,_i}, (121) 

and the corresponding space over which the smoother acts: 

5, = span [U{0P^}:^.,_.+r [Ji^' ^ ^p-^^l^-^ ■ 
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This set is used in the Borneniann-Erdmann-Kornhuber (BEK) refinement [111] 
and we utilize this set for the estimates of 3D local refinement. Since green refinement 
simply bisects a simplex, the modified basis function is the same as the one before 
the bisection due to linear interpolation. So the set of DOF in (121) corresponds 
to DOF created by red refinement and corresponding coarse DOF (father DOF). The 
following crucial result from [111] establishes a bound for the number of nodes used for 
smoothing. This indicates that the BPX preconditioner has provably optimal (linear) 
computational complexity per iteration on the resulting 3D mesh produced by the 
BEK refinement procedure. 

Lemma 4.1 The total number oj nodes used for smoothing satisfies the hound: 

Y.^i<\^J-\^^- (122) 

i=o 

Proof. See [111, Lemma 1]. □ 

The above lemma constitutes the first computational complexity optimality result 
in 3D for the BPX preconditioner as reported in [1<).']. A similar result for 2D red- 
green refinement was given by Oswald [11", page95]. In the general case of local 
smoothing operators which involve smoothing over newly created basis functions plus 
some additional set of local neighboring basis functions, one can extend the arguments 
from [111] and [110] using shape regularity. 



Hierarchical Basis Methods and Their Stabilizations 



HB methods exploit a 2-level hierarchical decomposition of the DOF. They are divided 
into the coarse (the ones inherited from previous levels) and the fine (the ones that 
are newly introduced) nodes. In fact, in the operator setting, this decomposition is a 
direct consequence of the direct decomposition of the finite element space as follows: 



^3 — ^j-i 



Hence, Aj can be represented by a two-by-two block form: 



where Aj 



A, 



^12 ! ^21 ' 



A 



aU) 

^12 



A. 



(123) 



and A22 correspond to coarse-coarse, coarse-fine, fine-coarse. 



and fine-fine discretization operators respectively. The same 2-level decomposition 
carries directly to the matrix setting. 

As mentioned earlier, HB methods suffer from the condition number growth. This 
makes makes HB methods especially ineffective in the 3D local refinement setting. 
As we mentioned earlier, wavelet-like stabilizations of the HB methods [lOI] provide 
a final solution to the condition number growth problem. The motivation for the 
stabilization hinges on the following idea. The BPX decomposition gives rise to basis 
functions which are not locally supported, but they decay rapidly outside a local 
support region. This allows for locally supported approximations as illustrated in 
Figures 3 and 4. 

The wavelet modified hierarchical basis (WMHB) methods [101, 112, 113] can 
be viewed as an approximation of the wavelet basis stemming from the BPX 
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Figure 3. Left: Hierarchical basis function without modification. Wavelet 
modified hierarchical basis functions. Middle: One iteration of symmetric Gauss- 
Seidel approximation. Right: One iteration of Jacobi approximation. 




Figure 4. Lower view of middle and left basis functions in Figure 3. 



decomposition [114]. A similar wavelet-like multilevel decomposition approach 
was taken in [115], where the orthogonal decomposition is formed by a discrete 
L2-equivalent inner product. This approach utilizes the same BPX two- level 
decomposition [UG, 115]. 

For local refinement setting, the other primary method of interest is the WMHB 
method. The WMHB methods can be described as additive or multiplicative Schwarz 
methods. In one of the previous papers [10-i], it was shown that the additive version 
of the WMHB method is optimal under certain types of red-green mesh refinement. 
Following the notational framework in [103, 105, 107, 113], this method is defined 
recursively as follows: 

Definition 4.1 The additive WMHB method D^^'' is defined for j — 1, . . . , J as 







B. 



U) 



with 
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With smooth PDE coefficients, optimal results were also established for the 
multiplicative version of the WMHB method in [f03, f07]. Our numerical experiments 
demonstrate such optimal results. This method can be written recursively as: 

Definition 4.2 The multiplicative WMHB method B'^^^ is defined as 







1) 



^12 



B. 



U) 



B. 





/ 



A 



'22 



a: 



(i) 



aU) 

^12 



B. 



with = 

i?22^ denotes an approximation of A22 , e.g. Gauss-Seidel or Jacobi approximation. 
For a more complete description of these and related algorithms, see [103, 107]. 



5. Practical Implementation of Fast Solvers 

The overall utility of any finite element code depends strongly on efficient 
implementation of its core algorithms and data structures. Finite element method 
becomes a viable tool in addressing realistic simulations only when these critical 
pieces come together. Theoretical results involving complexity are of little practical 
importance if the methods cannot be implemented. For algorithms involving 
data structures, this usually means striking a balance between storage costs and 
computational complexity. For instance, finding a minimal representation for a data 
set is only useful if the information can be accessed efficiently. We elaborate on the 
data structure and the implementation of the methods under consideration. 



5.1. Implementation of Hierarchical Basis Methods 



In HB methods, nodal basis functions are transformed into hierarchical basis functions 
via a nonsingular change of basis matrix: 



Y = 



I 

Y21 



Y12 
I + Y22 



where Y G Yr2 £ M^j-iX"^ Y21 G M"j><^^-i 

the number of fine DOF at level j as rij. Then, Nj ~ 



and Y22 e 
'^i + "-2 + 



We denote 
is the total 



number of DOF at level j. In this representation, we have assumed that the nodes 
are ordered with the nodes Nj-i inherited from the previous level listed first, and 
the rij new DOF listed second. For both wavelet modified (WMHB) and unmodified 
hierarchical basis (HB), the Y21 block represents the last Uj rows of the prolongation 
matrix; 

/ 



p3 



^21 



(124) 



In the HB case, the Y12 and 122 blocks are zero resulting in a very simple form: 

/ 

5^21 / 

Then, the original system (107) is related to the HB system through Y as follows. 

Ahb Uhb = Fhb, (125) 

= Y^AY, 
U - y C/hb, 
Fhh = Y^ F. 
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It will be shown later that the HBMG algorithm operates on sublocks of the HB 
system (125). We explicitly express the sublocks as follows. 







^ ^21 






A12 ' 




■ / " 
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A22 
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^2^1^22^21 



Ahb,12 — 
Ahb,21 — 
Ahb,22 = 



A 



12 



i21 



^2^1^22 
A22Y21 
A22, 



Next, we briefly describe the HBMG Algorithm. The HBMG routine can 
interpreted as an iterative process for solving the system (107) with an initial guess 
of Ui. 



Algorithm 3 



(i) smooth Ahb,22Uhb,2 = ^hb,2 - {Ahb,2iUhb,i) 
(a) form residual Rub.i = ^hh,i - {Ai^bAiUhb^) - ^hfc,i2C/h6,2 
(Hi) solve Ah6,iiJ7h6a = Rhb,i 

(iv) prolongate Uhb = Uhb + PUhb.i 

(v) smooth Ai^ba^Uhba = F^b,2 - (^h&.2iC^hfc,i) 

Smoothing involves the approximate solution of the linear system by a fixed number 
of iterations (typically one or two) with a method such as Gauss-Seidel, or Jacobi. In 
order to use HBMG as a preconditioner for CG, one has to make sure the pre-smoother 
is the adjoint of the post-smoother. One should also note that the algorithm can be 
simplified by first transforming the linear system (107) into the equivalent system 

A{U - Uj) ^F~AUj. 

In this setting, the initial guess is zero, and the HBMG algorithm recursively iterates 
towards the error with given residual on the right hand side. Simplification comes 
from the fact that terms in the parentheses are zero in Algorithm 3. 



5.1.1. The Gomputational and Storage Gomplexity of the HBMG method Utilizing 
the block structure of the stiffness matrix A^^^ in (123), the storage is in the following 
fashion: 
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JVj _ 1 X n j- 




^21 ^ ^ — > 




- aU) 
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rij xNj-i 






rij X rij 



where indicates a block stored rowwise, and x indicates a block which is not stored. 
By the symmetry in the bilinear form {DF{u)w,v) wrt w and v, ^ A21, hence 
A12 does not need to be stored. 
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Just like MG, the HBMG is an algebraic multilevel method as well. Namely, 
coarser stiffness matrices are formed algebraically through the use of variational 
conditions 

^ -'""^ .7 = 1,..., J; A:=A^''\ (126) 



where Pj_i denotes the prolongation operator from level j — 1 to j. Then, the only 

matrix to be stored for HB method is Pj-i, in which / is implicit, and therefore, does 
not have to be stored; 



pj 
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JVj_ixJVj_i 



Tlj X JVj _ 1 



In an adaptive scenario, new DOF are introduced in parts of the mesh where 

further correction or enrichment is needed. Naturally, the elements which are 
marked by the error estimator shrink in size by subdivision and the basis functions 
corresponding to fine nodes are forced to change more rapidly than the ones that 
correspond to the coarse nodes. Such rapidly changing components of the solution are 
described as the high frequency components. Smoothing is an operation which corrects 
the high frequency components of the solution, and is an integral part of the MG-like 
solvers. In MG, all DOF are exposed to smoothing which then requires to store all 
blocks of the stiffness matrix A. Unlike MG, the distinctive feature of the HBMG is 
that smoothing takes places only on basis functions corresponding to fine nodes. This 
feature leads us to make the crucial observation that V-cycle MG exactly becomes the 
HBMG method if smoothing is replaced by fine smoothing. One can describe HBMG 
style smoothing as fine smoothing. For fine smoothing, HBMG only needs to store the 
fine-fine interaction block A22 of A. It is exactly the fine smoothing that allows HB 
methods to have optimal storage complexity. 

In a typical case, nj is a small constant relative to Nj and has no relation to 
it. The HB method storage superiority stems from the fact that on every level j the 
storage cost is 0{nj). Fine smoothing is used for high frequency components, and 
this requires to store A22 block which is of size Uj x Uj. A22 is stored rowwise, and 
the storage cost is 0{nj). Coarse grid correction is used recursively for low frequency 
components, and this requires to store A12 and A21 blocks which are of size iVj_i x Uj 
and nj x Nj-\ respectively. Due to symmetry of the bilinear form, block is 
substituted by A21, and A21 is stored rowwise, and the storage cost is again 0{nj). 

Further strength of HB methods is that computational cost per cycle is 0{nj) 
on each level j. The preprocessing computational cost, namely computation of 
Ahb, 11, Ahb. 12, Ahb, 21, Ahb, 22 is 0{nj). Hence, in an adaptive refinement scenario, the 
overall computational complexity is achieved to be 0{N). 

Let us observe by a fictitious example how MG fails to maintain optimal storage 
complexity. Let us assume that the finest level is J, then N = Nj, Nj^i = N — nj, 
and count the total storage at each level. Here we take complexity constants to be 1 
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for simplicity. 

level J : TV 
level J - 1 : N - nj 

level 2 : N — nj — — ■ ■ ■ — 

level 1 : N — nj — nj-i — ■ ■ ■ ~ n-^ — n2. 

Adding up the cost at all levels, overall cost is: 

{J)N -{J- l)nj - (J - 2)nj_i 2n3 - 713. (127) 

Since nj's are small constants, then (127) < NJ, implying that overall storage is 
0{NJ). If numerous refinements are needed, or precisely, in an asymptotic scenario 
(i.e. J — > 00), the storage poses a severe problem. 



5.2. Sparse Matrix Products and the WMHB Implementation 

Our implementation relies on a total of four distinct sparse matrix data structures: 
compressed column (COL), compressed row (ROW), diagonal-row-column (DRC), 
and orthogonal-linked list (XLN). For detailed description of the data structures, the 
reader can refer to [1 0'l] . Each of these storage schemes attempts to record the location 
and value of the nonzeros using a minimal amount of information. The schemes differ 
in the exact representation which effects the speed and manner with which the data 
can be retrieved. XLN is an orthogonal-linked list data structure format which is the 
only dynamically "fillable" data structure used by our methods. By using variable 
length linked lists, rather than a fixed length array, it is suitable for situations where 
the total number of nonzeros is not known a priori. 

The key preprocessing step in the hierarchical basis methods, is converting the 
"nodal" matrices and vectors into the hierarchical basis. This operation involves sparse 
matrix- vector and matrix-matrix products for each level of refinement. To ensure that 
this entire operation has linear cost, with respect to the number of unknowns, the per- 
level change of basis operations must have a cost of 0{nj), where nj := Nj — iVj-i 
is the number of "new" nodes on level j. For the traditional multigrid algorithm this 
is not possible, since enforcing the variational conditions operates on all the nodes on 
each level, not just the newly introduced nodes. 

For WMHB, the Y12 and Y22 blocks are computed using the mass matrix, which 
results in the following formula: 

/ -inv [Mhb.ii] Mhb,i2 
Y21 I - r2iinv [Mhb.ii] A^hb,i2 

where the inv [•] is some approximation to the inverse which preserves the complexity. 
For example, it could be as simple as the inverse of the diagonal, or a low-order matrix 
polynomial approximation. The M^^ blocks are taken from the mass matrix in the 
HB basis: 

Mbb = YhbA^odaii^hb. (129) 

For the remainder of this section, we restrict our attention to the WMHB case. The 
HB case follows trivially with the two additional subblocks of Y set to zero. 



wmhb 



(128) 
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To reformulate the nodal matrix representation of the bilinear form in terms of 
the hierarchical basis, we must perform a triple matrix product of the form: 

^wmhb ~ ^ ^nodar 



In order to keep linear complexity, wc can only copy j4nodai fi fixed number of times, 
i.e. it cannot be copied on every level. Fixed size data structures are unsuitable for 
storing the product, since predicting the nonzero structure of ^^mhb i^^^ difficult 
as actually computing it. It is for these reasons that we have chosen the following 
strategy: First, copy ^nodai on the finest level, storing the result in an XLN which 
will eventually become Awmhb- Second, form the product pairwise, contributing the 
result to the XLN. Third, the last nj columns and rows of Awmhb are stripped off, 
stored in fixed size blocks, and the operation is repeated on the next level, using the 
All block as the new ^nodai^ 

Algorithm 4 (Wavelet Modified Hierarchical Change of Basis) 
Awmhb in XLN format. 



Copy Ai^Li 



While j >0 

(i) Multiply Awmhb = A^rnhhY as 



' An 


A12 ' 


+ = 


■ An 


A12 ■ 







yi2 " 


_ A21 


All . 


. A21 
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. ^21 


I22 



(ii) Multiply A^mbb = i^^Awmhb as 



An 
A21 



A12 
A22 



+ 



u -121 
-I12 -'2; 



An A12 
A21 A22 



(Hi) Remove A^i , A^^^ > ^22' blocks o/A^mbb storing in ROW, COL, and DRC formats 
respectively. 

(iv) After the removal, all that remains 0/ Awmhb is its A^{'^ block. 

(v) Let j = j - 1, descending to level 3 — 1. 

• End While. 

• Store the last Awmhb as Acoarse 

We should note that in order to preserve the complexity of the overall algorithm, 
all of the matrix-matrix algorithms must be carefully implemented. For example, the 
change of basis involves computing the products of An with Y12 and Y^^. To preserve 
storage complexity. Y12 must be kept in compressed column format, COL. For the 
actual product, the loop over the columns of F12 must be ordered first, then a loop 
over the nonzeros in each column, then a loop over the corresponding row or column 
in All. It iw cxac;tly for this reason, that one must be able to traverse An both by row 
and by column, which is why we have chosen an orthogonal-linked matrix structure 
for A during the change of basis (and hence An). 

To derive optimal complexity algorithms for the other products, it is enough 
to ensure that the outer loop is always over a dimension of size Uj. Due to the 
limited ways in which a sparse matrix can be traversed, the ordering of the remaining 
loops will usually be completely determined. Further gains can be obtained in the 
symmetric case, since only the upper or lower portion of the matrix needs to be 
explicitly computed and stored. 



(i) 
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6. The Finite Element ToolKit for the Einstein Constraints 

FETK [20] (see also [29, 30, 31]) is an adaptive multilevel finite element code in 
ANSI C developed by one of the authors (M.H.) over several years at Caltech and 
UC San Diego. It is designed to produce provably accurate numerical solutions to 
nonlinear covariant elliptic systems of tensor equations on 2- and 3-manifolds in an 
optimal or nearly-optimal way. FETK employs a posteriori error estimation, adaptive 
simplex subdivision, unstructured algebraic multilevel methods, global inexact Newton 
methods, and numerical continuation methods for the highly accurate numerical 
solution of nonlinear covariant elliptic systems on (Riemannian) 2- and 3-manifolds. 
In this section, we describe some of the key design features of FETK. A sequence of 
caxeful numerical examples producing initial data for an evolution simulation appear 
in [117]. 

6.1. The overall design o/FETK 

The finite element kernel library in FETK is referred to as MC (or Manifold Code). 
MC is an implementation of Algorithm 1, employing Algorithm 2 for nonlinear 
elliptic systems that arise in Step 1 of Algorithm 1. The linear Newton equations 
in each iteration of Algorithm 2 are solved with algebraic multilevel methods, and 
the algorithm is supplemented with a continuation technique when necessary. Several 
of the features of FETK are somewhat unusual, allowing for the treatment of very 
general nonlinear elliptic systems of tensor equations on domains with the structure 
of 2- and 3-manifolds. In particular, some of these features are: 

• Abstraction of the elliptic system: The elliptic system is defined only through 
a nonlinear weak form over the domain manifold, along with an associated 
linearization form, also defined everywhere on the domain manifold (precisely 
the forms {F{u),v) and {DF{u)w,v) in the discussions above). To use the a 
posteriori error estimator, a third function F{u) must also be provided (essentially 
the strong form of the problem) . 

• Abstraction of the domain manifold: The domain manifold is specified by giving a 
polyhedral representation of the topology, along with an abstract set of coordinate 
labels of the user's interpretation, possibly consisting of multiple charts. FETK 
works only with the topology of the domain, the connectivity of the polyhedral 
representation. The geometry of the domain manifold is provided only through 
the form definitions, which contain the manifold metric information, and through 
a oneChartO routine that the user provides to resolve chart boundaries. 

• Dimension independence: The same code paths are taken for both two- and three- 
dimensional problems (as well as for higher-dimensional problems). To achieve 
this dimension independence, FETK employs the simplex as its fundamental 
geometrical object for defining finite element bases. 

As a consequence of the abstract weak form approach to defining the problem, the 
complete definition of the constraints in the Einstein equations requires writing only 
1000 lines of C to define the two weak forms, and to define the oneChartO routine. 
Changing to a different problem, e.g., large deformation nonlinear elasticity, involves 
providing only a different definition of the forms and a different domain description. 
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6.2. Topology and geometry representation in FETK 

A datastmcture called the ringed-vertex (cf. [ ]) is used to represent meshes of d- 
simplices of arbitrary topology. This datastructure is illustrated in Figure 5. The 



Figure 5. Polyhedral manifold representation. The figure on the left shows two 
overlapping polyhedral (vertex) charts consisting of the two rings of simplices 
around two vertices sharing an edge. The region consisting of the two darkened 
triangles around the face / is denoted uj j , and represents the overlap of the two 
vertex charts. Polyhedral manifold topology is represented by FETK using the 
ringed vertex datastructure. The datastructure is illustrated for a given simplex 
s in the figure on the right; the topology primitives are vertices and d-simplices. 
The collection of the simplices which meet the simplex s at its vertices (which then 
includes those simplices that share faces as well) is denoted as lOg. (The set u)s 
includes s itself.) Edges are temporarily created during subdivision but are then 
destroyed (a similar ring datastructure is used to represent the edge topology). 



ringed- vertex datastructure is somewhat similar to the winged-edge, quad-edge, and 
edge-facet datastructures commonly used in the computational geometry community 
for representing 2-manifolds [118], but it can be used more generally to represent 
arbitrary d-manifolds, d = 2,3, .... It maintains a mesh of d-simplices with near 
minimal storage, yet for shape-regular (non-degenerate) meshes, it provides 0(l)-time 
access to all information necessary for refinement, un-refinement , and discretization 
of an elliptic operator. The ringed-vertex datastructure also allows for dimension 
independent implementations of mesh refinement and mesh manipulation, with one 
implementation covering arbitrary dimension d. An interesting feature of this 
datastructure is that the C structures used for vertices, simplices, and edges are all 
of fixed size, so that a fast array-based implementation is possible, as opposed to a 
less-efficient list-based approach commonly taken for finite element implementations 
on unstructured meshes. A detailed description of the ringed-vertex datastructure, 
along with a complexity analysis of various traversal algorithms, can be found in [ ]. 

Since FETK is based entirely on the d-simplex, for adaptive refinement it 
employs simplex bisection, using one of the simplex bisection strategies outlined 
earlier. Bisection is first used to refine an initial subset of the simplices in the mesh 
(selected according to some error estimates, discussed below), and then a closure 
algorithm is performed in which bisection is used recursively on any non-conforming 
simplices, until a conforming mesh is obtained. If it is necessary to improve element 
shape, FETK attempts to optimize the following simplex shape measure function for 
a given d-simplex s, in an iterative fashion, similar to the approach taken in [11!)]: 




ri{s,d) 




(130) 
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The quantity \s\ represents the (possibly negative) volume of the d-simplex, and \eij\ 
represents the length of the edge that connects vertex i to vertex j in the simplex. For 
d = 2 this is the shape-measure used in [119] with a slightly different normalization. 
For d = 3 this is the shape-measure developed in [120] again with a slightly different 
normalization. The shape measure above can be shown to be equivalent to the sphere 
ratio shape measure commonly used; the sphere measure is used in the black hole 
mesh generation algorithms we describe below. 

6.3. Discretization and adaptivity in FETK 

Given a nonlinear form {F(u),v), its linearization bilinear form {DF(u)w^v) , a 
Dirichlet function u, and collection of simplices representing the domain, FETK uses 
a default linear element to produce and then solve the implicitly defined nonlinear 
algebraic equations for the basis function coefficients in the expansion (87). The user 
can also provide his own element, specifying the number of degrees of freedom to 
be located on vertices, edges, faces, and in the interior of simplices, along with a 
quadrature rule, and the values of the basis functions at the quadrature points on the 
master element. Different element types may be used for different components of a 
coupled elliptic system. The availability of a user-defined general element makes it 
possible to, for example, use quadratic elements, which in the present setting would 
allow for the differentiation of the resulting solutions to the constraints for use as 
initial data in an evolution code. 

Once the equations are assembled and solved (discussed below), a posteriori error 
estimates are computed from the discrete solution to drive adaptive mesh refinement. 
The idea of adaptive error control in finite clement methods is to estimate the behavior 
of the actual solution to the problem using only a previously computed numerical 
solution, and then use the estimate to build an improved numerical solution by 
upping the polynomial order (p-refinement) or refining the mesh (^.-refinement) where 
appropriate. Note that this approach to adapting the mesh (or polynomial order) to 
the local solution behavior affects not only approximation quality, but also solution 
complexity: if a target solution accuracy can be obtained with fewer mesh points by 
their judicious placement in the domain, the cost of solving the discrete equations is 
reduced (sometimes dramatically) because the number of unknowns is reduced (again, 
sometimes dramatically). Generally speaking, if an elliptic equation has a solution 
with local singular behavior, such as would result from the presence of abrupt changes 
in the coefficients of the equation (e.g., the conformal metric components in the present 
case), then adaptive methods tend to give dramatic improvements over non-adaptive 
methods in terms of accuracy achieved for a given complexity price. Two examples 
illustrating bisection-based adaptivity patterns (driven by a completely geometrical 
"error" indicator simply for illustration) are shown in Figure 6. 

To drive the adaptivity in FETK, we employ a residual error estimator, based on 
computing an upper bound on the nonlinear residual as given in (110). Reference [2ii] 
contains the detailed derivation of the indicator (110) used in FETK for a general 
nonlinear elliptic system of tensor equations of the form (20)-(22). The error estimator 
provides a bound on the error in the VF^'^-norm, 1 < p < oo, for an approximation of 
the form (85) to the solution of the weak formulation (30). 
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Figure 6. Examples illustrating the 2D and 3D adaptive mesh refinement 
algorithms in FETK. The right-most figure in each row shows a close-up of the 
area where most of the refinement occurred in each example. 

6.4^. Solution of linear and nonlinear systems with FETK 

When a system of nonlinear finite element equations must be solved in FETK, the 
global inexact-Newton Algorithm 2 is employed, where the linearization systems 
are solved by linear multilevel methods. When necessary, the Newton procedure 
in Algorithm 2 is supplemented with a user-defined normalization equation for 
performing an augmented system continuation algorithm. The linear systems arising 
as the Newton equations in each iteration of Algorithm 2 are solved using a completely 
algebraic multilevel algorithm. Either refinement-generated prolongation matrices 
Pj-i, or user-defined prolongation matrices Pj^i in a standard YSMP-row-wise sparse 
matrix format, are used to define the multilevel hierarchy algebraically. In particular, 
once the single "fine" mesh is used to produce the discrete nonlinear problem F{u) = 
along with its linearization Au = f for use in the Newton iteration in Algorithm 2, a 
J-level hierarchy of linear problems is produced algebraically using the variational 
conditions (126) recursively. As a result, the underlying multilevel algorithm is 
provably convergent in the case of self-adjoint-positive matrices [121]. Moreover, the 
multilevel algorithm has provably optimal 0{N) convergence properties under the 
standard assumptions for uniform refinements [86], and is nearly-optimal 0(A^log A^) 
under very weak assumptions on adaptively refined problems [ ] . 

Coupled with the superlinear convergence properties of the outer inexact Newton 
iteration in Algorithm 2, this leads to an overall complexity of 0{N) or 0{N log N) for 
the solution of the discrete nonlinear problems in Step 1 of Algorithm 1. Combining 
this low-complexity solver with the judicious placement of unknowns only where 
needed due to the error estimation in Step 2 and the subdivision algorithm in Steps 3-6 
of Algorithm 1, leads to a very effective low-complexity approximation technique for 
solving a general class of nonlinear elliptic systems on 2- and 3-manifolds. 

6.5. Availability o/FETK 

A fully-functional MATLAB version of the MC kernel of FETK for domains with 
the structure of Riemannian 2-manifolds, called MCLite, is available under the GNU 
copy left license at the following website: 

http : //www . FEtk . ORG/ 
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MCLiTE employs the ringed- vertex datastructure and implements the same adaptivity 
and solution algorithms that are used in MC. A number of additional tools developed 
for use with MC and MCLite are also available under a GNU license at the site 
above, including MALOC (a Minimal Abstraction Layer for Object-oriented C) and 
SG (an OpenGL-based Xll/Win32 polygon display tool with a linear programming- 
based OpenGL-to-Postscript generator) . SG was used to generate the pictures of finite 
element meshes appearing in this paper. 

6.6. Tetrahedral Mesh Generation for Single or Binary Gompact Objects 

The general problem of generating a finite element mesh given a three dimensional 
domain is quite complex and has been the subject of intensive research over the past 
40 years (for a relatively old comprehensive overview the reader is referred to [ ]). 
However, for the purpose of generating a tetrahedral mesh on a domain with the 
geometry necessary to describe a binary collision between compact objects we do not 
need a completely general method. In this section we will describe a rather simple 
method for generating high quality coarse meshes suitable for the computation of the 
initial data describing binary systems. The only restriction of the method is that the 
geometry of the domain must have an axis of symmetry. We note that this restriction 
applies only to the domain geometry and not the physics (i.e., the source terms in the 
constraint equations do not have to have any symmetries). 




Figure 7. The planar triangulation step in the construction of a finite element 
mesh for a binary black hole collision. The figure shows the discretization of the 
entire domain; due to the mirror symmetry of the domain only the region a; > 
was meshed. 

A typical domain for a binary black hole collision consists of the interior of a 
large sphere (the outer boundary) with two smaller spheres removed (the throats of 
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the holes) . Since the exact placement of the outer sphere relative to the throats is not 
important we lose no generality in making the centers of the three spheres colinear. In 
this case the geometry is symmetric about the line joining the three centers. Choose 
this line to be the z axis. Then the intersection of the domain with the x-z plane 
consists of a circle with two smaller circles removed from its interior. The method 
consists of first meshing this planar domain with triangular elements and then rotating 
the elements around the z-axis to form triangular tori. These tori are then subdivided 
into tetrahedra in such a way as to produce a conforming mesh. 

We illustrate the method with an example. Let the outer boundary have radius 
100 and let the two throats be centered at z = ±4 with radii 2 and 1 respectively. 
The geometry is chosen for illustrative purposes. We discretize the outer circle of the 
planar domain with 20 line segments and the inner circles with 8 and 10 line segments 
respectively. Two views of the triangulation are shown in figures 7 and 8. 




Figure 8. A closcup of the mesh from Figure 7 showing the details of the mesh 
around the two holes. 

The triangulation algorithm one employs for this part of the mesh gen- 
eration is not important. In particular, one could use any one of a num- 
ber of off-the-shelf codes freely available on the internet (for example, see 
www-2 . OS . emu. edu/^quake/triangle .html). A list of commercial and free mesh 
generators can be found at 

www-users . inf ormatik. rwth-aachen. de/^roberts/ software .html. 
The triangulation code we have chosen to use is one developed by one of the authors 
(D.B.). We have the ability to produce highly graded meshes since the inner circles of 
the planar domain will in general be much small than the outer boundary circle. We 
would like to emphasize the distinguishing feature of the triangulation code. Although 
there are various codes that generate highly grades meshes with low quality, the code 
developed by D.B. is superior than those because it generates highly graded meshes 
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Figure 9. The second stage of mesh generation. One of the discretized tori is 
shown along with the planar mesh from figure 7. 

Once the planar triangulation is obtained the tetrahedra are constructed as 
follows. Each vertex of the triangulation defines a circle in the three-dimensional 
space by rotation about the z-axis. These circles are discretized by computing their 
circumference and dividing by a measure of the size of the desired tetrahedra. The 
measure we have used is the average distance to the neighbors of the vertex in the 
triangulation. For vertices with a neighbor on the axis of symmetry the discretization 
of the circle is limited to between five and seven line segments, regardless of the 
circumference or average neighbor distance. This limits the outdegree of the vertices 
which lie on the z-axis and prevents elements with unusually small angles from being 
formed there. Once each triangular torus is discretized along each of its edges it is a 
relatively simple matter to subdivide each torus into tetrahedra taking care that the 
elements thus formed are reasonably shaped and consistent with their neighbors. In 
figure 9 we show one of these tori and in figure 10 we show the surface of the final 
mesh. 

The quality of the meshes produced by this method is normally very good. One 
measure of the quality of a tetrahedron which is often used is the ratio of the radii of 
its circumscribing to inscribing spheres, a = R/3r. The normalization factor is such 
that a regular tetrahedron has a — 1 and it has been proven [120] that this is a global 
minimum, i.e., no tetrahedron has a < 1. The maximum aspect ratio in the mesh 
shown in figure 10 is 1.91 while the average is 1.26. By comparison the aspect ratio of 
the reference element with vertices located at (0,0,0), (1,0,0), (0,1,0), and (0,0,1) 
is 1.37. 

Finally, we note that although the example above describes the generation of a 
mesh for a binary black hole simulation the method is not restricted to this type of 
problem. We can also generate meshes suitable for any type of neutron star - black 
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Figure 10. The final mesh generated from the planar domain shown in figure 7. 
On the left the surface of the outer boundary is shown, on the right the surfaces 
of the two holes. 

hole collision by replacing one or both of the inner circles of the planar domain with a 
set of curves which are preserved on the mesh. The curves may represent (expected) 
isosurfaces of density or some other quantity. In Figures 11 and 12 we show the mesh 
generated for a collision between two neutron stars, one expected to be oblate and one 
expected to be prolate. Note that the initial data computed on such a mesh may or may 
not have stars of such shape, we have merely chosen a coarse mesh which we expect 
to reflect the nature of the solution to the constraints equations. The actual mesh 
which solves the numerical problem to a given tolerance may have a much different 
node distribution. 




Figure 11. A mesh suitable for the calculation of the collision between two stars, 
one initially prolate and one initially oblate. The figure shows shows the planar 
triangulation about the two stars. 
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Figure 12. The "surfaces" of tfie two stars represented by the mesh in Figure f 2 



6. 7. Computing Conformal Killing Vectors 

We describe a simple algorithm for computing conformal Killing vectors and their 
Petrov-Galerkin approximations. It can be shown that for Dirichlet and Robin 
boundary conditions the homogeneous version of the momentum constraint 

ba{LVY''=Q (131) 

has only the trivial solution — 0. However, using a pure Neumann condition 
removes the well-posedness of the problem and leads to a method for computing the 
conformal Killing vectors of 7ab. Suppose that K'^ satisfies (131). This implies 

= {ba{LKr\ V'')mM) = A{K\ y"), (132) 
V G £1 . This can be posed as the generalized eigenvalue problem 

Find K" e Hljj s. t. FiK", V) = eG{K'' , V), (133) 
\/ V £ HIjj, where 

F(if",F")= / 2^l{EK)ab{EVT''dx, (134) 

J M 

GiK^.V)^ [ DaK^bbV^dx. (135) 

J M 

If there exists a if" such that this holds for some e ^ then the operator F — eG 
arising from the form through the Ricsz representation theorem is singular. For a given 
domain with a tetrahedral decomposition, FETK can discretize (133) to produce the 
matrix versions of F and G which can then be given to a general eigenvalue package 
such as EISPACK. The result would be a set of eigenvalues and eigenvectors. The 
eigenvectors corresponding to the eigenvalue e = 2/3 form an (orthogonal) basis for 
the kernel of the discrete momentum operator, representing a discrete approximation 
to the space of conformal Killing vectors of ■^ab- 
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6.8. Gomputation of the ADM Mass on Adaptive Meshes 

Here we describe an alternative expression for the ADM mass that is more appropriate 
for adaptive methods than that usually used in numerical relativity. The ADM mass 
is normally computed by the covariant surface integral 

M= — - f haA^cj) ds (136) 
271- Jr 

where R indicates the integral is to be taken over the outer sphere only. The mass is 
actually defined as the limit i? — s- cx) and assumes that the equivalent energy associated 
with the conformal metric vanishes. However, if this is the case, and if the outer 
sphere is reasonably far away from the region where the source terms are varying 
rapidly then the adaptive nature of FETK will ensure that the outer sphere remains 
relatively coarsely meshed, even in cases where the solution is computed to a high 
degree of accuracy. This is because a 1/r variation in (j) can be computed accurately 
with only a few elements. Hence the outer boundary will remain coarsely meshed and 
the surface integral above will not be able to be estimated with any degree of accuracy. 
However, by converting (136) to a volume integral using the Gauss theorem, and by 
then evaluating it along with any resulting surface integrals over inner boundaries, it is 
possible to gain from the locally adapted mesh rather than lose from it. In particular, 
one has: 




(137) 



(a) (b) (c) 



Figure 13. An equatorial cut of (a) seven-block and (b) thirteen-block systems, 
(c) Grid dimensions for the seven-block system. 



Using the Hamiltonian constraint (14) this is 

M.-^/^(ii?0+l(tri.)V (138) 

-^CAab+{LW)abf^~' ~27Tpcl>-^^ dx (139) 

+ ^ / "aV> ds. (140) 

^■"^ JdM-R 

Normally the inner boundaries will be highly refined (although one can construct cases 
where this does not happen) and so this method will give much more accurate results 
than the simple use of (136). 
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Figure 14. (a) 1-d cut through the BriU wave conformal factor 4' fine for 
problem (a) along the x-axis. (b) Errors ^coarse - tpmedium, ^medium - ^fine, 
and pointwise convergence order on the 1-d cut along the x-axis. 



6.9. Brill waves initial data on multi-block domains 

As the final numerical experiment, we demonstrate the Brill wave initial data 
generated by using FETK on a multi-block spherical domain. The comprehensive 
treatment of this experiment can be found in [I IT]. 

In General Relativity, initial data on a spatial 3D-slice has to satisfy the 
Hamiltonian and momentum constraint equations [123], 

\7^{K'^ ~ g'^K) = 

where Kij and K are the extrinsic curvature of the 3D-slice and its trace, respectively, 
and ^R the Ricci scalar associated to the spatial metric gij. 

Brill waves [1-4] constitute a simple yet rich example of initial data in numerical 
relativity. In such case the extrinsic curvature of the slice is zero, and the above 
equations reduce to a single one, stating that the Ricci scalar has to vanish: 

^R = Q. (141) 

If the spatial metric is given up to one unknown function, Eq. (141) in principle 
allows to solve for such function and thus complete the construction of the initial 
data. The Brill equation is a special case of (141), where the 3-metric is expressed 
through the conformal transformation gij — ip'^gij of an unphysical metric gij, with 
an unknown conformal factor ip. Equation (141) then becomes [125]: 

+ l^)^ = (142) 

where R and V? are the Ricci scalar curvature and Laplacian of the unphysical metric 
gij, respectively. 

We work with two specific choices for g(p, z): 

(a) Holz' form [ ]: qnip, z) = aup^e"^ , with amplitude an — 0.5; 

(b) toroidal form: qt{p, z) — atp^ exp ^— 2°'' with amplitude at ~ 0.05, 

radius po — 5, width in p-direction ap = 3.0 and width in z-direction = 2.5. 




Here we will focus on the axisymmetric case with the conformal metric given in 
cylindrical coordinates by 

where q{p^ z) is a function satisfying the following conditions: 

(i) regularity at the axis: q{p = 0, z) = 0, f^|p=o — 0, 

(ii) asymptotic flatness: q{p, z)\r^^ < 0(l/r^), where r is the spherical radius 

The Hamiltonian constraint equation (141) becomes a second order elliptic PDE, 
which with asymptotically flat boundary conditions at r ^ oo takes the form 

~VV(p,^) + nP,2)V'(p,-)-0, (144) 

Af 

^Ir^c =1 + — + 0(l/r2), (145) 

2r 

with the potential V{p, z) given by 

V = - \{qpp + q'L)- 

We numerically solve this equation using FETK on the 13-patch multi-block 
spherical domain (see flgure 13). We use domain parameters Rout = 30, Rmed = 7, 
ttc = 1.5, and grid dimension ratios N : Nr^inner '■ Nr^outer = 2:3: 12. Our low- 
medium-high resolution triple is TV = 32, = 36 and iV = 40, except for pointwise 
convergence tests on the x-axis (see flgure 14), where we use = 16, = 24 and 
= 36 (since they all differ by powers of 1.5). 

7. Conclusion 



In this paper we have described the use of flnite element methods to construct 
numerical solutions to the initial value problem of general relativity. We flrst reviewed 
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the classical York conformal decomposition, and gave a basic framework for deriving 
weak formulations. We briefly outlined the notation used for the relevant function 
spaces, and gave a simple weak formulation example. We then derived an appropriate 
symmetric weak formulation of the coupled constraint equations, and summarized 
a number of basic theoretical properties of the constraints. We also derived the 
linearization bilinear form of the weak form for use with stability analysis or Newton- 
like numerical methods. A brief introduction to adaptive finite element methods for 
nonlinear elliptic systems was then presented, and residual-type error indicators were 
derived. We presented several general a priori error estimates from [2G, 27, 28] for 
general Galerkin approximations to solutions equations such as the momentum and 
Hamiltonian constraints. The numerical methods employed by MC were described 
in detail, including the finite element discretization, the residual-based a posteriori 
error estimator, the adaptive simplex bisection strategy, the algebraic multilevel 
solver, and the Newton-based continuation procedure for the solution of the nonlinear 
algebraic equations which arise. We described a mesh generation algorithm for 
modeling compact binary objects, outlined an algorithm for computing conformal 
Killing vectors, describes the numerical approximation of the ADM mass, and gave an 
example showing the use of MC for solution of the coupled constraints in the setting 
of a binary compact object collision. 

The implementation of these methods in the ANSI C finite element code named 
MC was discussed in detail, including descriptions of the algorithms and data 
structures it employs. MC was designed specifically for solving general second-order 
nonlinear elliptic systems of tensor equations on Riemannian manifolds with boundary. 
The key feature of MC which makes it particularly useful for relativity applications 
is the unusually high degree of abstraction with which it can be used. The user need 
only supply two functions (one for a linear problem) in the form of a short C code 
file. These functions are generally coded exactly as the weak form of the equation 
and its linearization are written down (our initial data constraint specification in MC 
is close to a cut-and-paste of equations (76) and (78) into a C source file) . The user 
does not have to provide the eUiptic system in discrete form as is usually required 
in finite difference implementations, and does not normally have to supply detailed 
coefficient information. In particular, the user provides only the two forms {F(u),v) 
and {DF{u)w,v). If a posteriori error estimation is to be used, then the user must 
provide a third function F{u), which is essentially the strong form of the differential 
equation as needed for the error estimator given in Section 3.4. MC is also able to 
handle systems on a manifold whose atlas has more than one chart. 
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